Fix: Added a fix in CLI and fragCounter.R to take reference when cram…

… files are used. Updated the NAMESPACE and DESCRIPTION with current paths

DESCRIPTION

@@ -5,11 +5,12 @@ Authors@R: person(given = "Marcin", family = "Imielinski", email = "mimielinski@
 Description: GC bias curve is determined by loess regression of read count by GC and mappability scores. Segmentation is done by circular binary segmentation (CBS) algorithm after getting tumor/normal ratios of corrected read counts.
 biocViews:
   bioc::3.10/GenomicRanges
+  bioc::3.10/GenomeInfoDb
 Remotes:
-  github::mskilab/gUtils,
-  github::mskilab/bamUtils,
-  github::mskilab/ffTrack,
-  github::mskilab/skidb
+  github::mskilab-org/gUtils,
+  github::mskilab-org/bamUtils,
+  github::mskilab-org/ffTrack,
+  github::mskilab-org/skidb
 Depends:
   R (>= 3.4)
 Imports:
@@ -31,4 +32,4 @@ Suggests:
 License: GPL-2
 Encoding: UTF-8
 LazyData: true
-RoxygenNote: 7.2.1.9000
+RoxygenNote: 7.2.3

NAMESPACE

@@ -9,6 +9,7 @@ export(correctcov_stub)
 export(fragCounter)
 export(make.blacklist)
 export(multicoco)
+import(GenomeInfoDb)
 import(GenomicRanges)
 import(gUtils)
 import(rtracklayer)

R/fragCounter.R

@@ -1,6 +1,7 @@
 #' @import GenomicRanges
 #' @import gUtils
 #' @import rtracklayer
+#' @import GenomeInfoDb
 #' @importFrom GenomeInfoDb seqlengths
 #' @importFrom data.table data.table fread rbindlist set setkey setkeyv setnames transpose as.data.table
 #' @importFrom stats cor loess predict quantile
@@ -303,7 +304,7 @@ fragCounter = function(bam, skeleton, cov = NULL, midpoint = TRUE, window = 200,
     cov$reads.corrected = coco(cov, mc.cores = 1, fields = c('gc', 'map'), iterative = T, exome = TRUE, imageroot = imageroot)$reads.corrected
 
   } else {
-    cov = PrepareCov(bam, cov = NULL, midpoint = midpoint, window = window, minmapq = minmapq, paired = paired, outdir, reference = reference)
+    cov = PrepareCov(bam, cov = NULL, reference = reference, midpoint = midpoint, window = window, minmapq = minmapq, paired = paired, outdir)
     cov = correctcov_stub(cov, gc.rds.dir = gc.rds.dir, map.rds.dir = map.rds.dir, chr.sub = chr.sub)
     cov$reads.corrected = multicoco(cov, numlevs = 1, base = max(10, 1e5/window), mc.cores = 1, fields = c('gc', 'map'), iterative = T, mono = T)$reads.corrected
   }
@@ -524,7 +525,8 @@ GC.fun = function(win.size = 200, twobitURL = '~/DB/UCSC/hg19.2bit', twobit.win
 #' @author Marcin Imielinski
 #' @param bam path to .bam file
 #' @param skeleton string Input data.table with intervals for which there is coverage data
-#' @param cov Path to existing coverage rds or bedgraph 
+#' @param cov Path to existing coverage rds or bedgraph
+#' @param reference reference file (recommended for CRAM, [NULL] 
 #' @param midpoint If TRUE only count midpoint if FALSE then count bin footprint of every fragment interval
 #' @param window window / bin size
 #' @param minmapq Minimal map quality
@@ -536,6 +538,7 @@ GC.fun = function(win.size = 200, twobitURL = '~/DB/UCSC/hg19.2bit', twobit.win
 #' @export
 
 PrepareCov = function(bam, skeleton, cov = NULL, reference = NULL, midpoint = TRUE, window = 200, minmapq = 1, paired = TRUE, outdir = NULL, exome = FALSE, use.skel = FALSE) {
+  library(GenomeInfoDb) # ADDED BY TANUBRATA: Forcefully loading GenomeInfoDb to Namespace since it is failing for cram files
   if (exome == TRUE){
 #    cov = bam.cov.exome(bam, chunksize = 1e6, min.mapq = 1)
     cov = bam.cov.skel(bam, skeleton, chunksize = 1e6, min.mapq = 1, use.skel = use.skel)
@@ -553,7 +556,7 @@ PrepareCov = function(bam, skeleton, cov = NULL, reference = NULL, midpoint = TR
         paired = TRUE
       }
       if (paired) {
-        cov = bamUtils::bam.cov.tile(bam, window = window, chunksize = 1e6, midpoint = TRUE, min.mapq = 1)  ## counts midpoints of fragments
+        cov = bamUtils::bam.cov.tile(bam, window = window, chunksize = 1e6, midpoint = TRUE, min.mapq = 1, reference = reference)  ## counts midpoints of fragments
       }
       else {
         file.type = bamorcram(bam)
@@ -953,7 +956,7 @@ bam.cov.skel = function(bam.file, skeleton, chunksize = 1e5, min.mapq = 1, verbo
 
   ref = ''
 
-  if (file.type == 'cramo')
+  if (file.type == 'cram')
   {
     sl = tryCatch(seqlengths(FaFile(reference)), error = function(e) NULL)
     if (is.null(sl))

inst/extdata/frag

@@ -28,7 +28,8 @@ if (!exists('opt'))
         make_option(c("-q", "--minmapq"), type = "integer", default = 1, help = "Minimal map quality"),
         make_option(c("-p", "--paired"), type = "logical", default = TRUE, help = "Is paired"),
         make_option(c("-e", "--exome"), type = "logical", default = FALSE, help = "Use exons as bins instead of fixed window"),
-        make_option(c("-r", "--chrsub"), type = "logical", default = TRUE, help = "remove chr prefix from seqnames"),
+	make_option(c("-r", "--reference"), type = "character", default = NULL, help = "Reference FASTA path (required for CRAM)"),
+        make_option(c("-a", "--chrsub"), type = "logical", default = TRUE, help = "remove chr prefix from seqnames"),
         make_option(c("-u", "--use.skel"), type = "logical", default = FALSE, help = "Use user defined regions instead of default exome skeleton"),
         make_option(c("-s", "--skeleton"), type = "character", help = "Path to skeleton file"),
         make_option(c("-o", "--outdir"), type = "character", default = './', help = "Directory to dump output into"),
@@ -51,8 +52,9 @@ if (!exists('opt'))
 suppressWarnings(suppressPackageStartupMessages(library(fragCounter)))
 suppressWarnings(suppressPackageStartupMessages(library(skidb)))
 suppressWarnings(suppressPackageStartupMessages(library(data.table)))
+suppressWarnings(suppressPackageStartupMessages(library(GenomeInfoDb)))
 message(dr.str)
 
-out = fragCounter(bam = opt$bam, cov = opt$cov, window = opt$window, gc.rds.dir = opt$gcmapdir, map.rds.dir = opt$gcmapdir, minmapq = opt$minmapq, paired = opt$paired, outdir = opt$outdir, exome = opt$exome, use.skel = opt$use.skel, skeleton = opt$skeleton, chr.sub = opt$chrsub)
+out = fragCounter(bam = opt$bam, cov = opt$cov, window = opt$window, gc.rds.dir = opt$gcmapdir, map.rds.dir = opt$gcmapdir, minmapq = opt$minmapq, paired = opt$paired, outdir = opt$outdir, exome = opt$exome, use.skel = opt$use.skel, skeleton = opt$skeleton, chr.sub = opt$chrsub, reference = opt$reference)
 
 saveRDS(out, paste(opt$outdir, 'cov.rds', sep = '/'))