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

 Type: Package

 Title: Copy number calling and SNV classification using

     targeted short read sequencing

-Version: 1.15.2

-Date: 2019-06-21

+Version: 1.15.3

+Date: 2019-08-10

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

         email="[email protected]"),

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

@@ -33,6 +33,7 @@ export(preprocessIntervals)

 export(processMultipleSamples)

 export(readCoverageFile)

 export(readCurationFile)

+export(readLogRatioFile)

 export(runAbsoluteCN)

 export(segmentationCBS)

 export(segmentationHclust)

@@ -5,6 +5,8 @@ NEW FEATURES

     o flag segments in poor quality regions

     o Dx.R now reports chromosomal instability scores

+    o Added readLogRatioFile function to read GATK4 DenoiseReadCounts

+      output files containing log2 tumor/normal ratios

 BUGFIXES

@@ -1099,3 +1099,20 @@ c(test.num.copy, round(opt.C))[i], prior.K, mapping.bias.ok, seg.id, min.variant

     return(c(ccf_point, ccf_lower, ccf_upper))

 }

+.getExonLrs <- function(seg.gr, tumor, log.ratio, idx = NULL) {

+    if (!is.null(idx)) {

+        tumor <- tumor[idx]

+        log.ratio <- log.ratio[idx]

+    }    

+    ov.se <- findOverlaps(seg.gr, tumor)

+    exon.lrs <- lapply(seq_len(length(seg.gr)), function(i) log.ratio[subjectHits(ov.se)[queryHits(ov.se) == 

+        i]])

+    exon.lrs <- lapply(exon.lrs, function(x) subset(x, !is.na(x) & !is.infinite(x)))

+    # estimate stand. dev. for target logR within targets. this will be used as proxy

+    # for sample error.

+    targetsPerSegment <- vapply(exon.lrs, length, integer(1))

+    if (!sum(targetsPerSegment > 50, na.rm = TRUE) && !is.null(idx)) {

+        .stopRuntimeError("Only tiny segments.")

+    }

+    exon.lrs

+}

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

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

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

-    x <- H5Fopen(file)

+    x <- H5Fopen(file, flags = "H5F_ACC_RDONLY")

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

new file mode 100755

@@ -0,0 +1,30 @@

+#' Read file containing interval-level log2 tumor/normal ratios 

+#' 

+#' Read log2 ratio file produced by external tools like The Genome Analysis 

+#' Toolkit version 4.

+#' 

+#' @param file Log2 coverage file.

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

+#' extension. Currently GATK4 DenoiseReadCounts format supported.

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

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

+#' @return A \code{GRange} with the log2 ratio.

+#' @author Markus Riester

+#' @examples

+#' 

+#' tumor.coverage.file <- system.file("extdata", "example_tumor.txt", 

+#'     package="PureCN")

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

+#' 

+#' @export readLogRatioFile

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

+    if (missing(format)) format <- "GATK4"

+    if (format == "GATK4") return(.readLogRatioFileGATK4(file, zero = FALSE))

+}

+

+.readLogRatioFileGATK4 <- function(file, zero = FALSE) {

+    x <- read.delim(file, comment.char = "@", as.is = TRUE)

+    gr <- GRanges(x[,1], IRanges(start = x[,2], end = x[,3]))

+    gr$log.ratio <- x[,4]

+    gr

+}

@@ -619,22 +619,20 @@ runAbsoluteCN <- function(normal.coverage.file = NULL,

     }

     # get target log-ratios for all segments

-    ov.se <- findOverlaps(seg.gr, tumor)

-    exon.lrs <- lapply(seq_len(nrow(seg)), function(i) log.ratio[subjectHits(ov.se)[queryHits(ov.se) == 

-        i]])

-    exon.lrs <- lapply(exon.lrs, function(x) subset(x, !is.na(x) & !is.infinite(x)))

-    

-    # estimate stand. dev. for target logR within targets. this will be used as proxy

-    # for sample error.

-    targetsPerSegment <- vapply(exon.lrs, length, integer(1))

-    

-    if (!sum(targetsPerSegment > 50, na.rm = TRUE)) {

-        .stopRuntimeError("Only tiny segments.")

-    }

+    exon.lrs <- .getExonLrs(seg.gr, tumor, log.ratio)

     sd.seg <- max(median(vapply(exon.lrs, sd, double(1)), na.rm = TRUE),

                   min.logr.sdev)

-   

+

+    if (sum(!tumor$on.target, na.rm = TRUE)) {

+        sd.seg.ontarget <- median(vapply(

+            .getExonLrs(seg.gr, tumor, log.ratio, idx = tumor$on.target), 

+            sd, double(1)), na.rm = TRUE)

+

+        sd.seg.offtarget <- median(vapply(

+            .getExonLrs(seg.gr, tumor, log.ratio, idx = !tumor$on.target), 

+            sd, double(1)), na.rm = TRUE)

+    } 

     # if user provided seg file, then we do not have access to the log-ratios and

     # need to use the user provided noise estimate also, don't do outlier smoothing

     # when we use already segmented data

@@ -669,6 +667,12 @@ runAbsoluteCN <- function(normal.coverage.file = NULL,

     if (sum(li < 0) > 0) 

         .stopRuntimeError("Some segments have negative size.")

     flog.info("Mean standard deviation of log-ratios: %.2f", sd.seg)

+    if (sum(!tumor$on.target, na.rm = TRUE)) {

+        flog.info("Mean standard deviation of on-target log-ratios only: %.2f",

+            sd.seg.ontarget)

+        flog.info("Mean standard deviation of off-target log-ratios only: %.2f",

+            sd.seg.offtarget)

+    }    

     log.ratio.offset <- rep(0, nrow(seg))

     flog.info("2D-grid search of purity and ploidy...")

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

         help = "Input: NormalDB.rds file. Generated by NormalDB.R."),

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

                 default = NULL, help = "Input: Segmentation file"),

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

+                default = NULL, help = "Input: Log2 copy number ratio file"),

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

         default = formals(PureCN::runAbsoluteCN)$sex[[2]],

         help = "Input: Sex of sample. ? (detect), diplod (non-diploid chromosomes removed), F or M [default %default]"),

@@ -136,6 +138,7 @@ if (!is.null(snp.blacklist)) {

 }

 seg.file <- opt$segfile

+log.ratio <- opt$logratiofile

 normalDB <- opt$normaldb

 sampleid <- opt$sampleid

 out <- opt[["out"]]

@@ -180,7 +183,6 @@ if (file.exists(file.rds) && !opt$force) {

     ret <- readCurationFile(file.rds)

     if (is.null(sampleid)) sampleid <- ret$input$sampleid

 } else {

-    tumor.coverage.file.orig <- tumor.coverage.file

     if (!is.null(normalDB)) {

         #if (!is.null(seg.file)) stop("normalDB and segfile do not work together.")

         normalDB <- readRDS(normalDB)

@@ -195,7 +197,8 @@ if (file.exists(file.rds) && !opt$force) {

                 normal.coverage.file <- calculateTangentNormal(tumor.coverage.file,

                     normalDB)

             }

-        } else if (is.null(normal.coverage.file) && is.null(seg.file)) {

+        } else if (is.null(normal.coverage.file) && is.null(seg.file) &&

+                   is.null(log.ratio)) {

             stop("Need either normalDB or normal.coverage.file")

         }

         normal.coverage.file

@@ -234,11 +237,20 @@ if (file.exists(file.rds) && !opt$force) {

         BPPARAM <- bpparam()

         flog.info("Using default BiocParallel backend. You can change the default in your ~/.Rprofile file.") 

     }

+    if (!is.null(log.ratio)) {

+        flog.info("Reading %s...", log.ratio)

+        log.ratio <- readLogRatioFile(log.ratio)

+        flog.info("Reading %s...", tumor.coverage.file)

+        tumor.coverage.file <- readCoverageFile(tumor.coverage.file)

+        tumor.coverage.file <- subsetByOverlaps(tumor.coverage.file, log.ratio)

+        log.ratio <- log.ratio$log.ratio

+    }    

     ret <- runAbsoluteCN(normal.coverage.file = normal.coverage.file,

             tumor.coverage.file = tumor.coverage.file, vcf.file = opt$vcf,

             sampleid = sampleid, plot.cnv = TRUE,

             interval.file = opt$intervals,

             genome = opt$genome, seg.file = seg.file,

+            log.ratio = log.ratio,

             test.purity = test.purity,

             test.num.copy = seq(0, opt$maxcopynumber),

             sex = opt$sex,

new file mode 100755

Binary files /dev/null and b/inst/extdata/example_gatk4_denoised_cr.tsv.gz differ

new file mode 100644

@@ -0,0 +1,34 @@

+% Generated by roxygen2: do not edit by hand

+% Please edit documentation in R/readLogRatioFile.R

+\name{readLogRatioFile}

+\alias{readLogRatioFile}

+\title{Read file containing interval-level log2 tumor/normal ratios}

+\usage{

+readLogRatioFile(file, format, zero = NULL)

+}

+\arguments{

+\item{file}{Log2 coverage file.}

+

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

+extension. Currently GATK4 DenoiseReadCounts format supported.}

+

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

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

+}

+\value{

+A \code{GRange} with the log2 ratio.

+}

+\description{

+Read log2 ratio file produced by external tools like The Genome Analysis 

+Toolkit version 4.

+}

+\examples{

+

+tumor.coverage.file <- system.file("extdata", "example_tumor.txt", 

+    package="PureCN")

+coverage <- readCoverageFile(tumor.coverage.file)

+

+}

+\author{

+Markus Riester

+}

new file mode 100644

@@ -0,0 +1,10 @@

+context("readLogRatioFile")

+

+test_that("Example data matches", {

+    logratio.file <- system.file("extdata", "example_gatk4_denoised_cr.tsv.gz", 

+        package = "PureCN")

+    logratio <- readLogRatioFile(logratio.file)

+    expect_equal(21, length(logratio))

+    expect_equal(0.109473, logratio$log.ratio[1], tolerance = .00001)

+    expect_equal(-0.185664, logratio$log.ratio[21], tolerance = .00001)

+})

@@ -552,6 +552,7 @@ Argument name      & Corresponding PureCN argument & PureCN function \\

 --normaldb         & normalDB (serialized with \Rfunction{saveRDS}) & 

     \Rfunction{calculateTangentNormal}, \Rfunction{filterTargets} \\

 --segfile          & seg.file           & \Rfunction{runAbsoluteCN}  \\

+--logratiofile     & log.ratio          & \Rfunction{runAbsoluteCN}  \\

 --sex              & sex                & \Rfunction{runAbsoluteCN}  \\

 --genome           & genome             & \Rfunction{runAbsoluteCN}  \\

 --intervals        & interval.file      & \Rfunction{runAbsoluteCN}   \\