diff --git a/R/getCNT.R b/R/getCNT.R index 7dab4b6..0b6a999 100644 --- a/R/getCNT.R +++ b/R/getCNT.R @@ -1,17 +1,17 @@ #' get the main copynumber state for each chromosome -#' +#' #' @param seq.dat the data.frame of sequencing data -#' +#' #' @details Get the main copynumber state for each chromosome. Calculate the total proportion of #' the chromosome associated with each unique copynumber state; return the copynumber state that is most #' represented. -#' +#' #' @return a data.frame containing the main copynumber state for each input chromosome -#' -#' @examples +#' +#' @examples #' seq.dat <- preprocessHRD( seq.dat ) #' CN.dat <- getCNt( seq.dat ) -#' +#' #' @export @@ -24,35 +24,44 @@ getCNt <- function( seq.dat ) # CN.dat (data.frame), then number { main.CN <- rep(0, length(unique(seq.dat$chromosome))) - + ct1 <- 0 for( chr in unique(seq.dat$chromosome)) { ct1 <- ct1 + 1 - + # reduce data to one chromosome at a time dat <- seq.dat[seq.dat$chromosome == chr,] - + # remove any NAs or this will crash dat <- dat[!is.na(dat$CNt),] CNt.frac <- rep(0, length(unique(dat$CNt))) ct2 <- 0 - + # get the fraction of the chromosome association with each copynumber state for( CNt in unique(dat$CNt) ) { ct2 <- ct2 + 1 CNt.frac[ct2] <- sum(dat$frac.chr[dat$CNt == CNt]) } - + # record values + associated CNt out <- data.frame(CNt=unique(dat$CNt), CNt.frac=CNt.frac) - + # get copynumber state most prominent in the chromosome - main.CN[ct1] <- out$CNt[which(CNt.frac == max(CNt.frac))] - + if(length(CNt.frac) > 0 && any(!is.na(CNt.frac))) { + max_frac_indices <- which(CNt.frac == max(CNt.frac, na.rm = TRUE)) + if(length(max_frac_indices) > 0) { + main.CN[ct1] <- out$CNt[max_frac_indices[1]] # Use the first if there are multiple maxima + } else { + main.CN[ct1] <- NA + } + } else { + main.CN[ct1] <- NA + } + } - + CN.dat <- data.frame(chromosome=unique(seq.dat$chromosome), main.CN=main.CN) return(CN.dat) } diff --git a/R/preprocessHRD.R b/R/preprocessHRD.R index 93ccd2a..a2d7c5e 100644 --- a/R/preprocessHRD.R +++ b/R/preprocessHRD.R @@ -1,12 +1,12 @@ #' Define additional properties of the sequencing data necessary to compute HRD -#' +#' #' @param seq.dat the raw sequencing data #' @param ref the reference genome to use (grch37, grch38) -#' @details define chromosome size, arms, telomere position, and centromere position from reference data; +#' @details define chromosome size, arms, telomere position, and centromere position from reference data; #' compute allelic imbalance and segment size -#' +#' #' @return seq.dat -#' +#' #'@examples #'seq.dat <- sub01.segments[ sub01.segments$chromosome == "chr1",] #' @@ -36,59 +36,65 @@ preprocessHRD <- function( seq.dat, ref ) print(paste(ref, "is not a valid reference genome.", sep = " ")) stop("select one of: grch37, grch38") } - + # check the sequencing data input to make sure it has the data we need if( any(!(seq.cols.needed %in% colnames(seq.dat)))) { - print(paste("column", + print(paste("column", seq.cols.needed[which(!(seq.cols.needed %in% colnames(seq.dat)))], "missing in the seq data.", sep=" ")) print(paste("Looking for columns:", seq.cols.needed, sep=" ")) print(paste("Found:", seq.cols.needed[seq.cols.needed %in% colnames(seq.dat)])) stop("Column mismatch. Exiting.") - } - - # remove rows w/ missing entries - if( any( is.na(seq.dat$CNt) )) - { - print(paste("Row: ", which(is.na(seq.dat$CNt)), " contains NA; removing.", sep = "")) - seq.dat <- subset(seq.dat, !is.na(seq.dat$CNt)) } - + + # fix chromosome numbers + seq.dat$chromosome[seq.dat$chromosome==23]<-"X" + seq.dat<-seq.dat[!seq.dat$chromosome %in% c(24,"Y"),] + # Check if any values in seq.dat$chromosome do NOT start with "chr" + if(any(!grepl("^chr", seq.dat$chromosome))) { + ref.dat$chromosome <- sub("^chr", "", ref.dat$chromosome) + } + levels(seq.dat$chromosome) <- levels(ref.dat$chromosome) - - seq.dat$frac.chr <- (seq.dat$end.pos - seq.dat$start.pos) / ref.dat$chr.size[match(seq.dat$chromosome, ref.dat$chromosome)] - + + if(nrow(seq.dat) > 0) { + seq.dat$frac.chr <- (seq.dat$end.pos - seq.dat$start.pos) / ref.dat$chr.size[match(seq.dat$chromosome, ref.dat$chromosome)] + } else { + print("No data available for calculating frac.chr.") + return(seq.dat) + } + # segment length seq.dat$brk.len <- 0 seq.dat <- combineSeg(seq.dat, 3e06) - + # matches reference data to the correct chromosome in the subject data key = match(seq.dat$chromosome, ref.dat$chromosome) - + seq.dat$seg.len <- seq.dat$end.pos - seq.dat$start.pos seq.dat$chr.size <- ref.dat$chr.size[match(seq.dat$chromosome, ref.dat$chromosome)] - + seq.dat$AI <- (seq.dat$A > seq.dat$B) & (seq.dat$A != 1) & (seq.dat$B != 0) - + seq.dat$start.arm <- seq.dat$end.arm <- rep("NA", dim(seq.dat)[1]) - + ind = seq.dat$start.pos - ref.dat$centromere.start[key] > 1 seq.dat$start.arm[ind] = "p" seq.dat$start.arm[!ind] = "q" - + ind = seq.dat$end.pos - ref.dat$centromere.end[key] > 1 seq.dat$end.arm[ind] = "p" seq.dat$end.arm[!ind] = "q" - + seq.dat$cross.arm <- seq.dat$start.arm != seq.dat$end.arm - - + + # define end telomeres seq.dat$post.telomere <- (seq.dat$start.pos - ref.dat$p.telomere.end[key] <= 1000) seq.dat$pre.telomere <- (ref.dat$q.telomere.start[key] - seq.dat$end.pos <= 1000) # --- - + return(seq.dat) } # ------------------------------------------------------------------------------------- #