@@ -2,8 +2,8 @@ Package: PureCN

 Type: Package

 Title: Copy number calling and SNV classification using

     targeted short read sequencing

-Version: 1.9.18

-Date: 2018-01-14

+Version: 1.9.19

+Date: 2018-01-15

 Authors@R: c(person("Markus", "Riester", role=c("aut", "cre"),

         email="[email protected]"),

     person("Angad P.", "Singh", role="aut"))

@@ -40,6 +40,7 @@ Imports:

     VGAM,

     edgeR,

     tools,

+    rhdf5,

     limma

 Suggests:

     PSCBS,

@@ -115,6 +115,7 @@ importFrom(graphics,strwidth)

 importFrom(graphics,symbols)

 importFrom(graphics,text)

 importFrom(gridExtra,grid.arrange)

+importFrom(rhdf5,H5Fopen)

 importFrom(rtracklayer,import)

 importFrom(stats,C)

 importFrom(stats,aggregate)

@@ -16,6 +16,7 @@

 #' @param index.file The bai index. This is expected without the .bai file

 #' suffix, see \code{?scanBam}.

 #' @param keep.duplicates Keep or remove duplicated reads.

+#' @param ... Additional parameters passed to \code{ScanBamParam}.

 #' @return Returns total and average coverage by intervals.

 #' @author Markus Riester

 #' @seealso \code{\link{preprocessIntervals}

@@ -37,7 +38,8 @@

 #' @importFrom Rsamtools headerTabix ScanBamParam scanBamFlag

 #'             scanBam scanFa scanFaIndex TabixFile

 calculateBamCoverageByInterval <- function(bam.file, interval.file,

-    output.file = NULL, index.file = bam.file, keep.duplicates = FALSE) {

+    output.file = NULL, index.file = bam.file, keep.duplicates = FALSE,

+    ...) {

     intervalGr <- readCoverageFile(interval.file)

     param <- ScanBamParam(what = c("pos", "qwidth", "flag"),

@@ -46,7 +48,8 @@ calculateBamCoverageByInterval <- function(bam.file, interval.file,

                              isNotPassingQualityControls = FALSE,

                              isSecondaryAlignment = FALSE,

                              isDuplicate = NA

-                     )

+                     ),

+                ...

              )

     xAll <- scanBam(bam.file, index = index.file, param = param)

@@ -5,10 +5,12 @@

 #' 

 #' @param file Target coverage file.

 #' @param format File format. If missing, derived from the file 

-#' extension. Currently only GATK DepthofCoverage and CNVkit formats 

-#' supported.

+#' extension. Currently GATK3 DepthofCoverage, GATK4 CollectFragmentCounts 

+#' (hdf5), and CNVkit formats supported.

 #' @param zero Start position is 0-based. Default is \code{FALSE}

 #' for GATK, \code{TRUE} for BED file based intervals.

+#' @param read.length For output formats which do not provide both counts 

+#' and total coverages, approximate them using the specified read length.

 #' @return A \code{data.frame} with the parsed coverage information.

 #' @author Markus Riester

 #' @seealso \code{\link{calculateBamCoverageByInterval}}

@@ -19,13 +21,16 @@

 #' coverage <- readCoverageFile(tumor.coverage.file)

 #' 

 #' @importFrom tools file_ext

+#' @importFrom rhdf5 H5Fopen

 #' @export readCoverageFile

-readCoverageFile <- function(file, format, zero=NULL) {

+readCoverageFile <- function(file, format, zero=NULL, read.length = 100) {

     if (missing(format)) format <- .getFormat(file)

     if (format %in% c("cnn", "cnr")) {

         targetCoverage <- .readCoverageCnn(file, zero, format)

+    } else if (format %in% c("hdf5")) {

+        targetCoverage <- .readCoverageGatk4(file, zero, format, read.length)

     } else {

-        targetCoverage <- .readCoverageGatk(file, zero)

+        targetCoverage <- .readCoverageGatk3(file, zero)

     }

     .checkLowCoverage(targetCoverage)

     .checkIntervals(targetCoverage)

@@ -34,10 +39,11 @@ readCoverageFile <- function(file, format, zero=NULL) {

 .getFormat <- function(file) {

     ext <- file_ext(file)

     if (ext %in% c("cnn", "cnr")) return(ext)

-    "GATK"    

-}    

+    if (ext %in% c("hdf5", "h5")) return("hdf5")

+    "GATK"

+}

-.readCoverageGatk <- function(file, zero) {

+.readCoverageGatk3 <- function(file, zero) {

     if (!is.null(zero)) flog.warn("zero ignored for GATK coverage files.")

     inputCoverage <- utils::read.table(file, header = TRUE)

     if (is.null(inputCoverage$total_coverage)) inputCoverage$total_coverage <- NA

@@ -55,6 +61,20 @@ readCoverageFile <- function(file, format, zero=NULL) {

     targetCoverage

 }

+.readCoverageGatk4 <- function(file, zero, format, read.length) {

+    if (!is.null(zero)) flog.warn("zero ignored for GATK coverage files.")

+    x <- H5Fopen(file)

+    intervals <- data.frame(x$intervals$transposed_index_start_end)

+    intervals[, 1] <- x$intervals$indexed_contig_names[intervals[, 1] + 1]

+    targetCoverage <- GRanges(intervals[, 1], IRanges(intervals[, 2], intervals[, 3]))

+    targetCoverage$counts <- x$counts$values[,1]

+    # TODO

+    targetCoverage$coverage <- targetCoverage$counts * read.length * 2

+    targetCoverage <- .addAverageCoverage(targetCoverage)

+    targetCoverage$on.target <- TRUE

+    targetCoverage

+}

+

 .addAverageCoverage <- function(x) {

     x$average.coverage <- x$coverage/width(x)

     x

@@ -20,6 +20,8 @@ option_list <- list(

     make_option(c("--keepduplicates"), action = "store_true", 

         default = formals(PureCN::calculateBamCoverageByInterval)$keep.duplicates, 

         help = "Count reads marked as duplicates [default %default]"),

+    make_option(c("--removemapq0"), action = "store_true", default = FALSE, 

+        help = "Not count reads marked with mapping quality 0 [default %default]"),

     make_option(c("--outdir"), action = "store", type = "character", 

         default = NULL,

         help = "Output directory to which results should be written"),

@@ -73,7 +75,7 @@ interval.file <- normalizePath(interval.file, mustWork = TRUE)

 }

 getCoverageBams <- function(bamFiles, indexFiles, outdir, interval.file, 

-    force = FALSE, cpu = 1, keep.duplicates = FALSE) {

+    force = FALSE, cpu = 1, keep.duplicates = FALSE, removemapq0 = FALSE) {

     bamFiles <- bamFiles

     indexFiles <- indexFiles

@@ -102,7 +104,8 @@ getCoverageBams <- function(bamFiles, indexFiles, outdir, interval.file,

         } else {

             PureCN::calculateBamCoverageByInterval(bam.file = bam.file,

                 interval.file = interval.file, output.file = output.file,

-                index.file = index.file, keep.duplicates = keep.duplicates)

+                index.file = index.file, keep.duplicates = keep.duplicates,

+                mapqFilter = if (removemapq0) 1 else NA)

         }

         output.file

     }

@@ -121,7 +124,8 @@ getCoverageBams <- function(bamFiles, indexFiles, outdir, interval.file,

 coverageFiles <- NULL

 indexFiles <- NULL

-flog.info("Loading PureCN...")

+flog.info("Loading PureCN %s...", Biobase::package.version("PureCN"))

+

 suppressPackageStartupMessages(library(PureCN))

 if (!is.null(bam.file)) {

@@ -142,7 +146,7 @@ if (!is.null(bam.file)) {

     }    

     coverageFiles <- getCoverageBams(bamFiles, indexFiles, outdir, 

-        interval.file, force, opt$cpu, opt$keepduplicates) 

+        interval.file, force, opt$cpu, opt$keepduplicates, opt$removemapq0) 

 }

 ### GC-normalize coverage -----------------------------------------------------

@@ -58,7 +58,7 @@ if (is.null(genome)) stop("Need --genome")

     file.path(outdir, paste0(prefix, assay, genome, suffix))

 }

-flog.info("Loading PureCN...")

+flog.info("Loading PureCN %s...", Biobase::package.version("PureCN"))

 if (length(coverageFiles)) {

     output.file <- .getFileName(outdir,"normalDB",".rds", assay, genome)

     if (file.exists(output.file) && !opt$force) {

@@ -149,7 +149,7 @@ if (!is.null(file.rds) && file.exists(file.rds)) {

 normalizePath(dirname(out), mustWork = TRUE)

-flog.info("Loading PureCN...")

+flog.info("Loading PureCN %s...", Biobase::package.version("PureCN"))

 suppressPackageStartupMessages(library(PureCN))

 library(futile.logger)

new file mode 100755

Binary files /dev/null and b/inst/extdata/example_normal5.hdf5 differ

@@ -5,7 +5,7 @@

 \title{Function to calculate coverage from BAM file}

 \usage{

 calculateBamCoverageByInterval(bam.file, interval.file, output.file = NULL,

-  index.file = bam.file, keep.duplicates = FALSE)

+  index.file = bam.file, keep.duplicates = FALSE, ...)

 }

 \arguments{

 \item{bam.file}{Filename of a BAM file.}

@@ -20,6 +20,8 @@ the \code{\link{readCoverageFile}} function.}

 suffix, see \code{?scanBam}.}

 \item{keep.duplicates}{Keep or remove duplicated reads.}

+

+\item{...}{Additional parameters passed to \code{ScanBamParam}.}

 }

 \value{

 Returns total and average coverage by intervals.

@@ -4,17 +4,20 @@

 \alias{readCoverageFile}

 \title{Read coverage file}

 \usage{

-readCoverageFile(file, format, zero = NULL)

+readCoverageFile(file, format, zero = NULL, read.length = 100)

 }

 \arguments{

 \item{file}{Target coverage file.}

 \item{format}{File format. If missing, derived from the file 

-extension. Currently only GATK DepthofCoverage and CNVkit formats 

-supported.}

+extension. Currently GATK3 DepthofCoverage, GATK4 CollectFragmentCounts 

+(hdf5), and CNVkit formats supported.}

 \item{zero}{Start position is 0-based. Default is \code{FALSE}

 for GATK, \code{TRUE} for BED file based intervals.}

+

+\item{read.length}{For output formats which do not provide both counts 

+and total coverages, approximate them using the specified read length.}

 }

 \value{

 A \code{data.frame} with the parsed coverage information.

@@ -53,3 +53,16 @@ test_that("CNVkit *cnr example data is parsed correctly", {

     expect_equal(coverage$on.target, c(FALSE, FALSE, FALSE, FALSE, 

         TRUE))

 })

+

+test_that("GATK4 *hdf5 example data is parsed correctly", {

+    coverageFile <- system.file("extdata", "example_normal5.hdf5", 

+        package = "PureCN")

+    coverage <- readCoverageFile(coverageFile)

+    expect_equal(length(coverage), 10)

+    expect_equal(head(start(coverage)), 

+        c(3598833, 3599562, 3607444, 3624039, 3638537, 3639872))

+    expect_equal(head(coverage$counts), 

+        c(127, 305, 78, 699, 566, 344))

+

+    expect_equal(head(coverage$on.target), rep(TRUE, 6))

+})