Merge branch 'refactor_inputs' into dev

mskilab-org · May 11, 2024 · f27a76c · f27a76c
2 parents 1a063cc + 79186bd
commit f27a76c
Show file tree

Hide file tree

Showing 24 changed files with 2,772 additions and 545 deletions.
diff --git a/bin/Events.R b/bin/Events.R
@@ -39,7 +39,7 @@ withAutoprint(
               min.span = 1e6,
               max.small = 1e4)
     newtic = tic
-    overwritefun("newtic", "tic", "gGnome")
+    # overwritefun("newtic", "tic", "gGnome")
 
 
     ## call complex events

diff --git a/bin/Fusions.R b/bin/Fusions.R
@@ -0,0 +1,50 @@
+withAutoprint(
+    {
+        library(optparse)
+        options(bitmapType = "cairo")
+
+        options(error = function() {
+            traceback(2)
+            quit("no", 1)
+        })
+        if (!exists("opt")) {
+            option_list <- list(
+                make_option(c("-i", "--id"), type = "character", help = "sample id"),
+                make_option(c("-g", "--gGraph"), type = "character", help = "an RDS file contains a gGraph or JaBbA graph with cn annotation on nodes and edges"),
+                make_option(c("-r", "--gencode"), type = "character", help = "an RDS or GTF file of GENCODE"),
+                make_option(c("-o", "--outdir"), type = "character", default = "./", help = "Directory to dump output into"),
+                make_option(c("--cores"), type = "integer", default = 1L, help = "Number of cores")
+            )
+            parseobj <- OptionParser(option_list = option_list)
+            opt <- parse_args(parseobj)
+            saveRDS(opt, paste(opt$outdir, "cmd.args.rds", sep = "/"))
+        }
+
+        library(gGnome)
+        library(gUtils)
+        library(parallel) ## needed for mc.cores
+
+        ## setDTthreads(10)
+        if (grepl(".rds$", opt$gencode)) {
+            gencode <- readRDS(as.character(opt$gencode))
+        } else {
+            gencode <- rtracklayer::import(opt$gencode)
+        }
+
+
+        ## call complex events
+        ## fus = fusions(gG(jab = opt$gGraph), gencode, verbose = TRUE, opt$cores)
+        fus <- fusions(gG(jab = opt$gGraph), gencode, verbose = TRUE, mc.cores = opt$cores)
+
+        ## update events with sample id
+        if (length(fus)) {
+            fus$set(id = opt$id)
+            fus$set(mincn = fus$eval(edge = min(cn, na.rm = TRUE)))
+        }
+
+        saveRDS(fus, paste0(opt$outdir, "/", "fusions.rds"))
+
+        quit("no", 0)
+    },
+    echo = FALSE
+)
diff --git a/bin/ascat_seg.R b/bin/ascat_seg.R
@@ -1,6 +1,6 @@
 {
     library(optparse)
-    
+
     ## DO NOT FAIL SILENTLY
     options(error = function() {traceback(2); quit("no", 1)})
 
@@ -35,7 +35,7 @@
         opt = parse_args(parseobj)
 
         print(opt)
-        
+
         print(.libPaths())
         options(error=function()traceback(2))
 
@@ -59,7 +59,7 @@
     #'
     #' @param gg (gGraph) diploid gGraph (expect cn.low and cn.high as node metadata)
     #' @param verbose (logical) default FALSE
-    #' 
+    #'
     #' @return gGraph with field allele and cn
     diploid2haploid = function(gg, verbose = FALSE) {
         og.nodes.gr = gg$nodes$gr[, c("cn.low", "cn.high", "var.low", "var.high", "cn", "node.id", "nhets")]
@@ -286,7 +286,7 @@
                 ##                                                          log.p = T) +
                 ##                                                  pnbinom(y, mu = centers[tot.cn - j + 1],
                 ##                                                          size = centers[tot.cn - j + 1],
-                ##                                                          log.p = T))) 
+                ##                                                          log.p = T)))
                 ll = sapply(0:(floor(tot.cn/2)), function(j) sum(ppois(x,centers[j+1], log.p = T) +
                                                                 ppois(y,centers[tot.cn-j+1],log.p = T)))
                 ll = ll - min(ll)
@@ -300,7 +300,7 @@
                         ##                                                          log.p = T) +
                         ##                                                  pnbinom(y, mu = centers[tot.cn - j],
                         ##                                                          size = centers[tot.cn - j] / 2,
-                        ##                                                          log.p = T))) 
+                        ##                                                          log.p = T)))
                         ll = sapply(0:(floor((tot.cn - 1)/2)), function(j) sum(ppois(x,centers[j+1], log.p = T) +
                                                                     ppois(y,centers[tot.cn-j],log.p = T)))
                         ll = ll - min(ll)
@@ -315,7 +315,7 @@
                     ##                                                          log.p = T) +
                     ##                                                  pnbinom(y, mu = centers[tot.cn - j + 2],
                     ##                                                          size = centers[tot.cn - j + 2],
-                    ##                                                          log.p = T))) 
+                    ##                                                          log.p = T)))
                     ll = sapply(0:(floor((tot.cn + 1)/2)), function(j) sum(ppois(x,centers[j+1], log.p = T) +
                                                                             ppois(y,centers[tot.cn-j+2],log.p = T)))
                     ll = ll - min(ll)
@@ -363,7 +363,7 @@
 
         ## ## ###########
         ## ## phasing
-        ## ## ########### 
+        ## ## ###########
 
         ## ## iterate through all reference junctions and apply (wishful thinking) heuristic
         ## ##
@@ -561,7 +561,7 @@
     #' @param agt.fname (character) path to sites.txt
     #' @param min.frac (numeric) between 0 and 1, min frequency in normal to count as het site
     #' @param max.frac (numeric) between 0 and 1, max frequency in normal to count as het site
-    #' 
+    #'
     #' @return allele gTrack
     grab.hets = function(agt.fname = NULL,
                         min.frac = 0.2,
@@ -605,7 +605,7 @@
         if (!file.exists(agt.fname)) {
             stop("invalid file")
         }
-        
+
         gr = readRDS(agt.fname)
 
         if (!inherits(gr, "GRanges")) {
@@ -748,6 +748,8 @@
                             median(ratio, na.rm = T)] > 0.8,
                     "XX",
                     "XY")
+    # handle single chr case as XY
+    gender = ifelse(!is.na(gender), gender, "XY")
     message("The gender of this sample: ", gender)
 
     message("Starting ASCAT!!")

diff --git a/bin/cbsFH.R b/bin/cbsFH.R
@@ -8506,7 +8506,7 @@ ppgrid = function(
   else ## only row, only go left right
       M = (NLLc < NLLcl & NLLc < NLLcr)[-c(1, nrow(NLLc)), -c(1, ncol(NLLc)), drop = FALSE]
 
-  if (length(M)>1)
+  if (length(M)>1 & any(M, na.rm=TRUE))
       {
           ix = which(M, arr.ind= T);
           if (nrow(ix)>1)