The changes I have made this far

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)
 }

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)
 }
 # ------------------------------------------------------------------------------------- #