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

 Type: Package

 Title: Copy number calling and SNV classification using

     targeted short read sequencing

-Version: 1.5.32

-Date: 2017-01-03

+Version: 1.5.33

+Date: 2017-01-06

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

         email="[email protected]"),

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

@@ -13,7 +13,7 @@ Description: This package estimates tumor purity, copy number, loss of

     with standard somatic variant detection pipelines, and has support for

     tumor samples without matching normal samples.

 Depends:

-    R (>= 3.2),

+    R (>= 3.3),

     DNAcopy,

     VariantAnnotation (>= 1.14.1)

 Imports:

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

     Biostrings,

     rtracklayer,

     ggplot2,

+    futile.logger,

     edgeR,

     limma

 Suggests:

@@ -55,13 +55,22 @@ importFrom(S4Vectors,queryHits)

 importFrom(S4Vectors,subjectHits)

 importFrom(SummarizedExperiment,rowRanges)

 importFrom(data.table,data.table)

+importFrom(futile.logger,appender.tee)

+importFrom(futile.logger,flog.appender)

+importFrom(futile.logger,flog.fatal)

+importFrom(futile.logger,flog.info)

+importFrom(futile.logger,flog.threshold)

+importFrom(futile.logger,flog.warn)

 importFrom(ggplot2,aes)

 importFrom(ggplot2,aes_string)

 importFrom(ggplot2,element_text)

 importFrom(ggplot2,facet_wrap)

+importFrom(ggplot2,geom_boxplot)

+importFrom(ggplot2,geom_hline)

 importFrom(ggplot2,geom_line)

 importFrom(ggplot2,geom_point)

 importFrom(ggplot2,ggplot)

+importFrom(ggplot2,labs)

 importFrom(ggplot2,scale_alpha_continuous)

 importFrom(ggplot2,stat_density2d)

 importFrom(ggplot2,theme)

@@ -6,6 +6,7 @@ Changes in version 1.6.0

 - Polished plots, added new GC-normalization and volcano plots

 - Support for cell lines

 - Improved somatic vs. germline status calling

+- Contamination rate estimation.

 - Better mapping bias estimation and correction

 - Better copy number normalization using multiple best normals

 - New GC-normalization for smaller gene panels

@@ -16,6 +17,8 @@ Changes in version 1.6.0

 - Faster post.optimize=TRUE by not optimizing poor fits.

 - Tweaks to segmentationPSCBS

 - Lots of improvements to command line scripts.

+- Code cleanups (switch from inlinedocs to roxygen, from message/warn to

+  futile.logger)

 API CHANGES

@@ -40,6 +43,14 @@ API CHANGES

         normalized so that w[1] is 1

     - Removed ... from runAbsoluteCN

     - added centromeres to segmentation function.

+    - contamination.cutoff now specifies fraction of SNPs necessary to

+      call sample contaminated. Does not remove anymore, since these

+      are obviously informative for contamination rate estimation in the

+      PureCN model.

+    - Removed verbose from most functions, since messages are now controlled

+      with futile.logger.  

+    - Smoothing of log-ratios now optionally done by runAbsoluteCN, not

+      segmentation function.  

 PLANNED FEATURES 

@@ -50,3 +61,13 @@ PLANNED FEATURES

   MYC)

 - Switch to S4 data structures (maybe)

 - Whole dataset visualizations (maybe PureCN 1.8)

+

+

+ROADMAP

+

+- January 15: API freeze (no changes to existing API)

+- February 1: Defaults freeze (no tuning of default parameters anymore)

+- February 15: Feature freeze (only bugfixes, code cleanups and speedups

+  and documentation improvements) 

+- March 15: Code freeze (only bugfixes and documentation improvements)

+- April: Release

@@ -64,7 +64,7 @@ c(test.num.copy, round(opt.C))[i], prior.K))

 }    

 .calcSNVLLik <- function(vcf, tumor.id.in.vcf, ov, p, test.num.copy, 

-    C.posterior, C, opt.C, snv.model, prior.somatic, mapping.bias, snv.lr, 

+    C.likelihood, C, opt.C, snv.model, prior.somatic, mapping.bias, snv.lr, 

     sampleid = NULL, cont.rate = 0.01, prior.K, max.coverage.vcf, non.clonal.M,

     model.homozygous=FALSE, error=0.001, max.mapping.bias=0.8) {

@@ -84,9 +84,9 @@ c(test.num.copy, round(opt.C))[i], prior.K))

         haploid.penalty <- 2

     }

-    subclonal <- apply(C.posterior[queryHits(ov), ], 1, which.max) == ncol(C.posterior)

+    subclonal <- apply(C.likelihood[queryHits(ov), ], 1, which.max) == ncol(C.likelihood)

-    seg.idx <- which(seq_len(nrow(C.posterior)) %in% queryHits(ov))

+    seg.idx <- which(seq_len(nrow(C.likelihood)) %in% queryHits(ov))

     sd.ar <- sd(unlist(geno(vcf)$FA[, tumor.id.in.vcf]))

     # Fit variants in all segments

@@ -104,13 +104,13 @@ c(test.num.copy, round(opt.C))[i], prior.K))

         shape1 <- ar_all * dp_all + 1

         shape2 <- (1 - ar_all) * dp_all + 1

-        lapply(seq(ncol(C.posterior)), function(k) {

+        lapply(seq(ncol(C.likelihood)), function(k) {

             Ci <- c(test.num.copy, opt.C[i])[k]       

             priorM <- log(Ci + 1 + haploid.penalty)

             #priorM <- log(abs(2-Ci) + 1)

             priorHom <- ifelse(model.homozygous, -log(3), log(0))

-            skip <- test.num.copy > Ci | C.posterior[i, k] <=0

+            skip <- test.num.copy > Ci | C.likelihood[i, k] <=0

             p.ar <- lapply(c(0, 1), function(g) {

                 cns <- test.num.copy

@@ -166,8 +166,8 @@ c(test.num.copy, round(opt.C))[i], prior.K))

     likelihoods <- do.call(rbind, 

         lapply(seq_along(xx), function(i) Reduce("+", 

-            lapply(seq(ncol(C.posterior)), function(j) 

-                exp(xx[[i]][[j]]) * C.posterior[seg.idx[i], j]))))

+            lapply(seq(ncol(C.likelihood)), function(j) 

+                exp(xx[[i]][[j]]) * C.likelihood[seg.idx[i], j]))))

     colnames(likelihoods) <- c(paste("SOMATIC.M", test.num.copy, sep = ""), 

         paste("GERMLINE.M", test.num.copy, sep = ""), "GERMLINE.CONTHIGH", 

@@ -418,7 +418,7 @@ c(test.num.copy, round(opt.C))[i], prior.K))

 .getGoF <- function(result) {

     if (is.null(result$SNV.posterior$beta.model)) return(0)

     r <- result$SNV.posterior$beta.model$posteriors

-    e <- (r$ML.AR-r$AR)^2

+    e <- (r$ML.AR-r$AR.ADJUSTED)^2

     maxDist <- 0.2^2

     r2 <- max(1-mean(e,na.rm=TRUE)/maxDist,0)

     return(r2)

@@ -583,7 +583,7 @@ c(test.num.copy, round(opt.C))[i], prior.K))

 .optimizeGrid <- function(test.purity, min.ploidy, max.ploidy, test.num.copy = 0:7, 

-    exon.lrs, seg, sd.seg, li, max.exon.ratio, max.non.clonal, verbose = FALSE, debug = FALSE) {

+    exon.lrs, seg, sd.seg, li, max.exon.ratio, max.non.clonal) {

     ploidy.grid <- seq(min.ploidy, max.ploidy, by = 0.2)

     if (min.ploidy < 1.8 && max.ploidy > 2.2) {

         ploidy.grid <- c(seq(min.ploidy, 1.8, by = 0.2), 1.9, 2, 2.1, seq(2.2, max.ploidy, 

@@ -592,16 +592,12 @@ c(test.num.copy, round(opt.C))[i], prior.K))

     mm <- lapply(test.purity, function(p) {

         b <- 2 * (1 - p)

         log.ratio.offset <- 0

-        if (debug) 

-            message(paste(b, log.ratio.offset))

         lapply(ploidy.grid, function(D) {

             dt <- p/D

             llik.all <- lapply(seq_along(exon.lrs), function(i) .calcLlikSegmentExonLrs(exon.lrs[[i]], 

                 log.ratio.offset, max.exon.ratio, sd.seg, dt, b, D, test.num.copy))

             subclonal <- vapply(llik.all, which.max, double(1)) == 1

             subclonal.f <- length(unlist(exon.lrs[subclonal]))/length(unlist(exon.lrs))

-            if (debug) 

-                message(paste(sum(subclonal), subclonal.f))

             if (subclonal.f > max.non.clonal) 

                 return(-Inf)

             llik <- sum(vapply(llik.all, max, double(1)))

@@ -857,6 +853,7 @@ c(test.num.copy, round(opt.C))[i], prior.K))

     msg <- paste0(msg, "\n\nThis is most likely a user error due to",

         " invalid input data or parameters (PureCN ", 

         packageVersion("PureCN"), ").")

+    flog.fatal(paste(strwrap(msg),"\n"))

     stop( paste(strwrap(msg),"\n"), call.= FALSE)

 }

 .stopRuntimeError <- function(...) {

@@ -864,25 +861,41 @@ c(test.num.copy, round(opt.C))[i], prior.K))

     msg <- paste0(msg, "\n\nThis runtime error might be caused by",

         " invalid input data or parameters. Please report bug (PureCN ", 

         packageVersion("PureCN"), ").")

+    flog.fatal(paste(strwrap(msg),"\n"))

     stop( paste(strwrap(msg),"\n"), call.= FALSE)

 }

-

+.logHeader <- function(l) {

+    flog.info(strrep("-", 60))

+    flog.info("PureCN %s", as.character(packageVersion("PureCN")))

+    flog.info(strrep("-", 60))

+    idxSupported <- sapply(l, function(x) class(eval(x))) %in% 

+        c("character", "numeric", "NULL", "list", "logical") &

+        sapply(l, function(x) length(unlist(eval(x)))) < 12

+

+    l <- c(l[idxSupported],lapply(l[!idxSupported], function(x) "<data>"))

+    argsStrings <- paste(sapply(seq_along(l), function(i) 

+        paste0("-", names(l)[i], " ", paste(eval(l[[i]]),collapse=","))),

+        collapse=" ")

+    flog.info("Arguments: %s", argsStrings)

+}

+.logFooter <- function() {

+    flog.info("Done.")

+    flog.info(strrep("-", 60))

+}    

 .calcGCmetric <- function(gc.data, coverage) { 

     gcbins <- split(coverage$average.coverage, gc.data$gc_bias < 0.5); 

     mean(gcbins[[1]], na.rm=TRUE)/mean(gcbins[[2]], na.rm=TRUE) 

 }

-.checkGCBias <- function(normal, tumor, gc.data, max.dropout, verbose) {

+.checkGCBias <- function(normal, tumor, gc.data, max.dropout) {

     gcMetricNormal <- .calcGCmetric(gc.data, normal)

     gcMetricTumor <- .calcGCmetric(gc.data, tumor)

-    if (verbose) { 

-        message("AT/GC dropout: ", round(gcMetricTumor, digits=2),

-        " (tumor), ", round(gcMetricNormal, digits=2), " (normal).")

-    }    

+    flog.info("AT/GC dropout: %.2f (tumor), %.2f (normal). ", 

+        gcMetricTumor, gcMetricNormal)

     if (gcMetricNormal < max.dropout[1] || 

         gcMetricNormal > max.dropout[2] ||

         gcMetricTumor  < max.dropout[1] ||

         gcMetricTumor  > max.dropout[2]) {

-        warning("High GC-bias in normal or tumor. Is data GC-normalized?")

+        flog.warn("High GC-bias in normal or tumor. Is data GC-normalized?")

         return(TRUE)

     }

     return(FALSE)

@@ -10,7 +10,6 @@

 #' function.

 #' @param tumor Tumor coverage read in by the \code{\link{readCoverageGatk}}

 #' function.

-#' @param verbose Verbose output.

 #' @return \code{numeric(nrow(tumor))}, tumor vs. normal copy number log-ratios

 #' for all targets.

 #' @author Markus Riester

@@ -22,21 +21,17 @@

 #'     package="PureCN")

 #' normal <- readCoverageGatk(normal.coverage.file)

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

-#' log.ratio <- calculateLogRatio(normal, tumor, verbose=FALSE)

+#' log.ratio <- calculateLogRatio(normal, tumor)

 #' 

 #' @export calculateLogRatio

-calculateLogRatio <- function(normal, tumor, verbose=TRUE) {

+calculateLogRatio <- function(normal, tumor) {

     # make sure that normal and tumor align

     if (!identical(as.character(normal[, 1]), as.character(tumor[, 1]))) {

         .stopUserError("Interval files in normal and tumor different.")

     }

-    if (verbose) {

-        message("Average coverage: ", 

-            round(mean(tumor$average.coverage, na.rm=TRUE), digits=0), 

-            "X (tumor), ",

-            round(mean(normal$average.coverage, na.rm=TRUE), digits=0),

-            "X (normal).")

-    }    

+    flog.info("Average coverage: %.0fX (tumor) %.0fX (normal).", 

+        mean(tumor$average.coverage, na.rm=TRUE), 

+        mean(normal$average.coverage, na.rm=TRUE))

     total.cov.normal <- sum(as.numeric(normal$coverage), na.rm = TRUE)

     total.cov.tumor <- sum(as.numeric(tumor$coverage), na.rm = TRUE)

@@ -57,10 +57,10 @@ log.ratio.cutoffs = c(-0.9, 0.9), failed = NULL, all.genes = FALSE) {

     bm <- res$results[[id]]$SNV.posterior$beta.model

     if (!is.null(bm)) {

-        segids <- bm$posterior$seg.id

+        segids <- bm$posteriors$seg.id

         calls$num.snps.segment <- sapply(calls$seg.id, function(i) 

             sum(segids==i,na.rm=TRUE))

-        calls$loh <- bm$posterior$ML.M.SEGMENT[match(calls$seg.id, segids)] == 0 

+        calls$loh <- bm$posteriors$ML.M.SEGMENT[match(calls$seg.id, segids)] == 0 

     }

     calls <- calls[, !grepl("^\\.",colnames(calls))]

@@ -1,5 +1,5 @@

 # Make CMD check happy

-globalVariables(names=c("gcIndex", "gcNum", "..level.."))

+globalVariables(names=c("..level.."))

 #' Correct for GC bias

 #' 

@@ -110,7 +110,7 @@ plot.max.density = 50000) {

         if (density == "Low") {

             print(ggplot(gcPlot, aes_string(x="gc_bias", y="average.coverage")) + 

                 geom_point(color='red', alpha=0.2) + 

-                geom_line(data = plotMed, aes(x = gcIndex, y = gcNum), color = 'blue') + 

+                geom_line(data = plotMed, aes_string(x = 'gcIndex', y = 'gcNum'), color = 'blue') + 

                 xlab("GC content") + ylab("Coverage") + 

                 theme(axis.text = element_text(size= 6), axis.title = element_text(size=16)) + 

                 facet_wrap(~ norm_status, nrow=1))

@@ -119,7 +119,7 @@ plot.max.density = 50000) {

             geom_point(color="blue", alpha = 0.1) + 

             stat_density2d(aes(fill = ..level..), geom="polygon") + 

             scale_alpha_continuous(limits=c(0.1, 0), breaks=seq(0, 0.1, by = 0.025)) + 

-            geom_line(data = plotMed, aes(x = gcIndex,y = gcNum), color = 'red') + 

+            geom_line(data = plotMed, aes_string(x = 'gcIndex',y = 'gcNum'), color = 'red') + 

             xlab("GC content") + ylab("Coverage") + 

             theme(axis.text = element_text(size = 16), axis.title = element_text(size = 16)) + 

             facet_wrap(~norm_status, nrow=1))

@@ -46,7 +46,7 @@ max.mean.coverage = 100, ... ) {

     normals.m <- scale(normals.m, 1/z, center=FALSE)

     normals.pca <- prcomp(t(normals.m[idx,]), ...)

-    sex.determined <- sapply(normals,getSexFromCoverage, verbose=is.null(sex))

+    sex.determined <- sapply(normals,getSexFromCoverage)

     if (is.null(sex)) {

         sex <- sex.determined

     } else {   

@@ -59,9 +59,8 @@ max.mean.coverage = 100, ... ) {

         for (i in seq_along(sex.determined)) {

             if (!is.na(sex.determined[i]) && sex[i] != "diploid" &&

                 sex.determined[i] != sex[i]) {

-                warning("Sex mismatch in ", normal.coverage.files[i], 

-                    ". Sex provided is ", sex, ", but could be ", 

-                    sex.determined[i])

+                flog.warn("Sex mismatch in %s. Sex provided is %s, but could be %s.", 

+                    normal.coverage.files[i], sex[i], sex.determined[i])

             }    

         }    

     }

@@ -109,7 +108,6 @@ createExonWeightFile <- function() {

 #' (>20) to estimate target log-ratio standard deviations. Should not overlap

 #' with files in \code{tumor.coverage.files}.

 #' @param target.weight.file Output filename.

-#' @param verbose Verbose output.

 #' @return A \code{data.frame} with target weights.

 #' @author Markus Riester

 #' @examples

@@ -127,12 +125,12 @@ createExonWeightFile <- function() {

 #' 

 #' @export createTargetWeights

 createTargetWeights <- function(tumor.coverage.files, normal.coverage.files,

-target.weight.file, verbose = TRUE) {

-    if (verbose) message("Loading coverage data...")

+target.weight.file) {

+    flog.info("Loading coverage data...")

     tumor.coverage <- lapply(tumor.coverage.files,  readCoverageGatk)

     normal.coverage <- lapply(normal.coverage.files,  readCoverageGatk)

     lrs <- lapply(tumor.coverage, function(tc) sapply(normal.coverage, 

-            function(nc) calculateLogRatio(nc, tc, verbose=verbose)))

+            function(nc) calculateLogRatio(nc, tc)))

     lrs <- do.call(cbind, lrs)

@@ -23,7 +23,6 @@

 #' \code{\link{createNormalDatabase}}.

 #' @param normalDB.min.coverage Exclude targets with coverage lower than 20

 #' percent of the chromosome median in the pool of normals.

-#' @param verbose Verbose output.

 #' @return \code{logical(length(log.ratio))} specifying which targets should be

 #' used in segmentation.

 #' @author Markus Riester

@@ -55,7 +54,7 @@

 #' @export filterTargets

 filterTargets <- function(log.ratio, tumor, gc.data, seg.file, 

     filter.lowhigh.gc = 0.001, min.targeted.base = 5, normalDB = NULL,

-    normalDB.min.coverage = 0.2, verbose) {

+    normalDB.min.coverage = 0.2) {

     # NA's in log.ratio confuse the CBS function

     targetsUsed <- !is.na(log.ratio) & !is.infinite(log.ratio) 

     # With segmentation file, ignore all filters

@@ -64,12 +63,12 @@ filterTargets <- function(log.ratio, tumor, gc.data, seg.file,

     if (!is.null(gc.data)) {

         .checkFraction(filter.lowhigh.gc, "filter.lowhigh.gc")

         targetsUsed <- .filterTargetsLowHighGC(targetsUsed, tumor,

-            gc.data, filter.lowhigh.gc, verbose)

+            gc.data, filter.lowhigh.gc)

     }

     targetsUsed <- .filterTargetsNormalDB(targetsUsed, tumor, normalDB,

-        normalDB.min.coverage, verbose)

+        normalDB.min.coverage)

     targetsUsed <- .filterTargetsTargetedBase(targetsUsed, tumor,

-        min.targeted.base, verbose)

+        min.targeted.base)

 }

 .checkNormalDB <- function(tumor, normalDB) {

@@ -90,7 +89,7 @@ filterTargets <- function(log.ratio, tumor, gc.data, seg.file,

 }

 .filterTargetsNormalDB <- function(targetsUsed, tumor, normalDB,

-normalDB.min.coverage, verbose) {

+normalDB.min.coverage) {

     if (is.null(normalDB)) return(targetsUsed)

     # make sure that normalDB matches tumor

     if (!.checkNormalDB(tumor, normalDB)) {

@@ -113,40 +112,40 @@ normalDB.min.coverage, verbose) {

     nAfter <- sum(targetsUsed)

-    if (verbose && nAfter < nBefore) { 

-        message("Removing ", nBefore-nAfter, " targets with low coverage ",

-            "in normalDB.")

+    if (nAfter < nBefore) { 

+        flog.info("Removing %i targets with low coverage in normalDB.", 

+            nBefore-nAfter)

     }

     targetsUsed

 }

-.filterTargetsChrHash <- function(targetsUsed, tumor, chr.hash, verbose) {

+.filterTargetsChrHash <- function(targetsUsed, tumor, chr.hash) {

     if (is.null(chr.hash)) return(targetsUsed)

     nBefore <- sum(targetsUsed)

     targetsUsed <- targetsUsed & tumor$chr %in% chr.hash$chr

     nAfter <- sum(targetsUsed)

-    if (verbose && nAfter < nBefore) { 

-        message("Removing ", nBefore-nAfter, " targets on chromosomes ",

-            "outside chr.hash.")

+    if ( nAfter < nBefore) { 

+        flog.info("Removing %i targets on chromosomes outside chr.hash.", 

+            nBefore-nAfter)

     }

     targetsUsed

 }

-.filterTargetsTargetedBase <- function(targetsUsed, tumor, min.targeted.base,

-    verbose) {

+.filterTargetsTargetedBase <- function(targetsUsed, tumor, min.targeted.base) {

     if (is.null(min.targeted.base)) return(targetsUsed)

     nBefore <- sum(targetsUsed)

     targetsUsed <- targetsUsed & !is.na(tumor$targeted.base) & 

         tumor$targeted.base >= min.targeted.base

     nAfter <- sum(targetsUsed)

-    if (verbose && nAfter < nBefore) message("Removing ", nBefore-nAfter, 

-        " small targets")

+    if (nAfter < nBefore) {

+        flog.info("Removing %i small targets.", nBefore-nAfter)

+    }    

     targetsUsed

 }

 .filterTargetsLowHighGC <- function(targetsUsed, tumor, gc.data,

-    filter.lowhigh.gc, verbose) {

+    filter.lowhigh.gc) {

     gc.data <- gc.data[match(as.character(tumor[,1]), gc.data[,1]),]

     qq <- quantile(gc.data$gc_bias, p=c(filter.lowhigh.gc, 

         1-filter.lowhigh.gc), na.rm=TRUE)

@@ -156,8 +155,8 @@ normalDB.min.coverage, verbose) {

         !(gc.data$gc_bias < qq[1] | gc.data$gc_bias > qq[2])

     nAfter <- sum(targetsUsed)

-    if (verbose && nAfter < nBefore) message("Removing ", 

-        nBefore-nAfter, " low/high GC targets")

-

+    if (nAfter < nBefore) {

+        flog.info("Removing %i low/high GC targets.", nBefore-nAfter)

+    }

     targetsUsed

 }

@@ -20,8 +20,9 @@

 #' contamination. If a matched normal is available, this value is ignored,

 #' because homozygosity can be confirmed in the normal.

 #' @param contamination.cutoff Count SNPs in dbSNP with allelic fraction

-#' smaller than the first value, if found on most chromosomes, remove all with

-#' AF smaller than the second value.

+#' smaller than the first value or greater than 1-first value, if found on most chromosomes, mark sample

+#' as contaminated if the fraction of putative contamination SNPs exceeds

+#' the second value.

 #' @param min.coverage Minimum coverage in tumor. Variants with lower coverage

 #' are ignored.

 #' @param min.base.quality Minimim base quality in tumor. Requires a \code{BQ}

@@ -40,7 +41,6 @@

 #' SNPs. Ignored in case a matched normal is provided in the VCF.

 #' @param interval.padding Include variants in the interval flanking regions of

 #' the specified size in bp. Requires \code{target.granges}.

-#' @param verbose Verbose output.

 #' @return A list with elements \item{vcf}{The filtered \code{CollapsedVCF}

 #' object.} \item{flag}{A flag (\code{logical(1)}) if problems were

 #' identified.} \item{flag_comment}{A comment describing the flagging.}

@@ -61,10 +61,10 @@

 #' @importFrom stats pbeta

 filterVcfBasic <- function(vcf, tumor.id.in.vcf = NULL, 

 use.somatic.status = TRUE, snp.blacklist = NULL, af.range = c(0.03, 0.97),

-contamination.cutoff = c(0.05,0.075), min.coverage = 15, min.base.quality = 25,

+contamination.cutoff = c(0.075,0.02), min.coverage = 15, min.base.quality = 25,

 min.supporting.reads = NULL, error = 0.001, target.granges = NULL,

-remove.off.target.snvs = TRUE, model.homozygous = FALSE, interval.padding = 50,

-verbose=TRUE) {

+remove.off.target.snvs = TRUE, model.homozygous = FALSE, 

+interval.padding = 50) {

     flag <- NA

     flag_comment <- NA

@@ -74,25 +74,73 @@ verbose=TRUE) {

     if (use.somatic.status) {

         n <- nrow(vcf)

         vcf <- .testGermline(vcf, tumor.id.in.vcf)

-        if (verbose) message(paste("Removing", n-nrow(vcf), 

-            "non heterozygous (in matched normal) germline SNPs."))

+        flog.info("Removing %i non heterozygous (in matched normal) germline SNPs.", 

+            n-nrow(vcf))

     } else {

         info(vcf)$SOMATIC <- NULL

     }    

+    # find supporting read cutoffs based on coverage and sequencing error

+    #--------------------------------------------------------------------------

     if (is.null(min.supporting.reads)) {

-        min.supporting.reads <- calculatePowerDetectSomatic(

-            mean(geno(vcf)$DP[,tumor.id.in.vcf], na.rm=TRUE),

-            purity=1,ploidy=2, error=error, verbose=FALSE)$k

-    }    

+        depths <- geno(vcf)$DP[,tumor.id.in.vcf]

+        minDepth <- round(log2(mean(depths)))

+

+        depths <- sort(unique(round(log2(depths))))

+        depths <- depths[depths>=minDepth]

+        depths <- 2^depths

+        cutoffs <- sapply(depths, function(d) 

+            calculatePowerDetectSomatic(d, ploidy=2, purity=1, error=error, 

+                verbose=FALSE)$k)

+        depths <- c(0,depths)

+    } else { 

+        depths <- 0

+        cutoffs <- min.supporting.reads

+    }

     n <- nrow(vcf)

-    # remove variants with insufficient reads. This includes reference alleles

-    # to remove homozygous germline variants.

-    vcf <- vcf[do.call(rbind, 

-        geno(vcf)$AD[,tumor.id.in.vcf])[,2] >= min.supporting.reads]

-    

+    .sufficientReads <- function(vcf, ref, depths, cutoffs) {

+        idx <- rep(TRUE, nrow(vcf))

+

+        .filterVcfByAD <- function(vcf, min.supporting.reads, depth) {

+            # remove variants with insufficient reads. This includes reference alleles

+            # to remove homozygous germline variants.

+            do.call(rbind, geno(vcf)$AD[,tumor.id.in.vcf])[,ifelse(ref,1,2)] >= 

+                min.supporting.reads | 

+                geno(vcf)$DP[,tumor.id.in.vcf] < depth

+           }

+

+        for (i in seq_along(cutoffs)) {

+            idx <- idx & .filterVcfByAD(vcf, cutoffs[i], depths[i])

+        }

+        idx

+    }    

+    idxNotReference <- .sufficientReads(vcf,ref=FALSE, depths, cutoffs)

+    vcf <- vcf[idxNotReference]

+    idxNotHomozygous <- .sufficientReads(vcf,ref=TRUE, depths, cutoffs)

+    #--------------------------------------------------------------------------

+

+    idx <- info(vcf)$DB & idxNotHomozygous &

+        ( unlist(geno(vcf)$FA[,tumor.id.in.vcf]) <  contamination.cutoff[1] | 

+        unlist(geno(vcf)$FA[,tumor.id.in.vcf]) > (1 - contamination.cutoff[1]))

+

+    fractionContaminated <- sum(idx)/sum(info(vcf)$DB)

+    if (fractionContaminated > 0) {

+    #    cntContPerChr <- sapply(seqlevelsInUse(vcf), function(seq) nrow(keepSeqlevels(vcf[idx],seq)))

+    #    cntAllPerChr <- sapply(seqlevelsInUse(vcf), function(seq) nrow(keepSeqlevels(vcf,seq)))

+    #    flog.info("Contamination p-value across chromosomes: %.2f.",

+    #        fisher.test(cntContPerChr, cntAllPerChr)$p.value)

+    }

+

+    # do we have many low allelic fraction calls that are in dbSNP on basically

+    # all chromosomes? then we found some contamination

+    if (sum(runLength(seqnames(rowRanges(vcf[idx])))>3) >= 20 &&

+        fractionContaminated >= contamination.cutoff[2]) {

+        flag <- TRUE

+        flag_comment <- "POTENTIAL SAMPLE CONTAMINATION"    

+    }

+

     # If we have a matched normal, we can distinguish LOH from homozygous

     # in 100% pure samples. If not we need to see a sufficient number

     # of alt alleles to believe a heterozygous normal genotype.    

@@ -100,20 +148,15 @@ verbose=TRUE) {

         af.range[2] <- 1

         if (model.homozygous) af.range[2] <- Inf

     } else {

-        vcf <- vcf[do.call(rbind, 

-            geno(vcf)$AD[,tumor.id.in.vcf])[,1] >= min.supporting.reads]

+        vcf <- vcf[idxNotHomozygous]

     }

     vcf <- vcf[unlist(geno(vcf)$DP[,tumor.id.in.vcf]) >= min.coverage]

     vcf <- vcf[unlist(geno(vcf)$FA[,tumor.id.in.vcf]) >= af.range[1]]

     # remove homozygous germline

     vcf <- vcf[!info(vcf)$DB | geno(vcf)$FA[,tumor.id.in.vcf] < af.range[2]]

-    if (verbose) message("Removing ", n-nrow(vcf), 

-        " SNPs with AF < ", af.range[1],

-        " or AF >= ", af.range[2], 

-        " or less than ", min.supporting.reads, 

-        " supporting reads or depth < ",

-        min.coverage, ".")

+    flog.info("Removing %i SNPs with AF < %.3f or AF >= %.3f or less than %i supporting reads or depth < %i.", 

+        n-nrow(vcf), af.range[1], af.range[2], cutoffs[1], min.coverage)

     if (!is.null(snp.blacklist)) {

         for (i in seq_along(snp.blacklist)) {

@@ -124,7 +167,7 @@ verbose=TRUE) {

             }

             n <- nrow(vcf)

             if (sum( rownames(vcf) %in% snp.blacklist.data[,1]) > 1  ) {

-                warning("Old SNP blacklists are deprecated. ", 

+                flog.warn("Old SNP blacklists are deprecated. %s", 

                     "Use either a BED file or a normal.panel.vcf.file.")

                 vcf <- vcf[!rownames(vcf) %in% snp.blacklist.data[,1],]

             } else {

@@ -132,79 +175,37 @@ verbose=TRUE) {

                 ov <- suppressWarnings(overlapsAny(vcf, blackBed))

                 vcf <- vcf[!ov]

             }    

-            if (verbose) message("Removing ", n-nrow(vcf), 

-                " blacklisted SNPs.")

-        }    

-    }

-    idx <- info(vcf)$DB & unlist(geno(vcf)$FA[,tumor.id.in.vcf]) < 

-        contamination.cutoff[1]

-

-    # do we have many low allelic fraction calls that are in dbSNP on basically

-    # all chromosomes? then we found some contamination

-    if (sum(runLength(seqnames(rowRanges(vcf[idx])))>1) >= 20) {

-        idx <- info(vcf)$DB & unlist(geno(vcf)$FA[,tumor.id.in.vcf]) < 

-            contamination.cutoff[2]

-        vcf <- vcf[which(!idx)]

-        if (verbose) message("Removing ", sum(idx, na.rm=TRUE), 

-            " contamination SNPs.")

-        flag <- TRUE

-        flag_comment <- "POTENTIAL SAMPLE CONTAMINATION"    

-    }

-    

-    # if we have a matched normal, we can filter homozygous germline SNPs

-    # easily, if not, we still should remove most of them by removing everything

-    # with AF >= 0.97. But sometimes samples are very noisy, for example due to

-    # contamination, and homozygous germline SNPs are between 0.9 and 1.  This

-    # code finds samples which have a lot of dbSNPs between 0.95 and 0.97 on

-    # almost all chromosomes.  If found, then remove everything that's in dbSNP

-    # above 0.9.

-    if (!use.somatic.status && !model.homozygous) {    

-        idx <- info(vcf)$DB & unlist(geno(vcf)$FA[,tumor.id.in.vcf]) > 

-            (1 - contamination.cutoff[1])

-

-        # do we have many high allelic fraction calls that are in dbSNP on

-        # basically all chromosomes?  then we found a very noisy sample

-        if (sum(runLength(seqnames(rowRanges(vcf[idx])))>2) >= 20) {

-            idx <- info(vcf)$DB & unlist(geno(vcf)$FA[,tumor.id.in.vcf]) > 

-                (1 - contamination.cutoff[2])

-            vcf <- vcf[which(!idx)]

-            if (verbose) message("Removing ", sum(idx, na.rm=TRUE), 

-                " noisy homozygous germline SNPs.")

-            flag <- TRUE

-            flag_comment <- .appendComment(flag_comment, 

-                "NOISY HOMOZYGOUS GERMLINE CALLS")    

+            flog.info("Removing %i blacklisted SNPs.", n-nrow(vcf))

         }    

     }

     if (!is.null(geno(vcf)$BQ)) {

-       n.vcf.before.filter <- nrow(vcf)

-       vcf <- vcf[which(as.numeric(geno(vcf)$BQ[,tumor.id.in.vcf])>=min.base.quality)]

-       if (verbose) message("Removing ", n.vcf.before.filter - nrow(vcf), 

-           " low quality variants with BQ<",min.base.quality,".") 

+        n.vcf.before.filter <- nrow(vcf)

+        vcf <- vcf[which(as.numeric(geno(vcf)$BQ[,tumor.id.in.vcf])>=min.base.quality)]

+        flog.info("Removing %i low quality variants with BQ < %i.", 

+            n.vcf.before.filter - nrow(vcf), min.base.quality) 

     } else if (!is.null(rowRanges(vcf)$QUAL)) {

-       n.vcf.before.filter <- nrow(vcf)

-       idx <- which(as.numeric(rowRanges(vcf)$QUAL)>=min.base.quality)

-       if (length(idx)) vcf <- vcf[idx]

-       if (verbose) message("Removing ", n.vcf.before.filter - nrow(vcf), 

-           " low quality variants with BQ<",min.base.quality,".") 

+        n.vcf.before.filter <- nrow(vcf)

+        idx <- which(as.numeric(rowRanges(vcf)$QUAL)>=min.base.quality)

+        if (length(idx)) vcf <- vcf[idx]

+        flog.info("Removing %i low quality variants with BQ < %i.", 

+            n.vcf.before.filter - nrow(vcf), min.base.quality) 

     }     

     if (!is.null(target.granges)) {

-        vcf <- .annotateVcfTarget(vcf, target.granges, interval.padding, verbose)

+        vcf <- .annotateVcfTarget(vcf, target.granges, interval.padding)

         if (remove.off.target.snvs) {

             n.vcf.before.filter <- nrow(vcf)

             # make sure all SNVs are in covered exons

             vcf <- vcf[info(vcf)$OnTarget>0]

-            if (verbose) message("Removing ", n.vcf.before.filter - nrow(vcf), 

-                " variants outside intervals.", 

-                " Set remove.off.target.snvs=FALSE to include.")

+            flog.info("Removing %i variants outside intervals.", 

+                n.vcf.before.filter - nrow(vcf))

         }        

     }

-    ##value<< A list with elements

     list(

-        vcf=vcf, ##<< The filtered \code{CollapsedVCF} object.

-        flag=flag, ##<< A flag (\code{logical(1)}) if problems were identified.

-        flag_comment=flag_comment ##<< A comment describing the flagging. 

+        vcf=vcf, 

+        flag=flag, 

+        flag_comment=flag_comment 

     )

 }

@@ -222,7 +223,6 @@ verbose=TRUE) {

 #' @param tumor.id.in.vcf The tumor id in the VCF file, optional.

 #' @param stats.file MuTect stats file.

 #' @param ignore MuTect flags that mark variants for exclusion.

-#' @param verbose Verbose output.

 #' @param \dots Additional arguments passed to \code{\link{filterVcfBasic}}.

 #' @return A list with elements \code{vcf}, \code{flag} and

 #' \code{flag_comment}.  \code{vcf} contains the filtered \code{CollapsedVCF},

@@ -244,15 +244,15 @@ filterVcfMuTect <- function(vcf, tumor.id.in.vcf = NULL, stats.file = NULL,

 ignore=c("clustered_read_position", "fstar_tumor_lod", "nearby_gap_events", 

 "poor_mapping_region_alternate_allele_mapq", "poor_mapping_region_mapq0", 

 "possible_contamination", "strand_artifact", "seen_in_panel_of_normals"),

-verbose=TRUE, ... ){

+... ){

     if (is.null(stats.file)) return(

-        filterVcfBasic(vcf, tumor.id.in.vcf, verbose=verbose, ...))

+        filterVcfBasic(vcf, tumor.id.in.vcf, ...))

     stats <- read.delim(stats.file, as.is=TRUE, skip=1)

     if (is.null(stats$contig) || is.null(stats$position)) {

         warning("MuTect stats file lacks contig and position columns.")

-        return(filterVcfBasic(vcf, tumor.id.in.vcf, verbose=verbose, ...))

+        return(filterVcfBasic(vcf, tumor.id.in.vcf, ...))

     }    

     gr.stats <- GRanges(seqnames=stats$contig, 

@@ -265,13 +265,13 @@ verbose=TRUE, ... ){

         n <- nrow(vcf)

         stats <- stats[subjectHits(ov),]

         vcf <- vcf[queryHits(ov)]

-        warning("MuTect stats file and VCF file do not align perfectly. ",

-         "Will remove ", n-nrow(vcf), " unmatched variants.")

+        flog.warn("MuTect stats file and VCF file do not align perfectly. Will remove %i unmatched variants.",

+            n-nrow(vcf))

     }    

     if (is.null(stats$failure_reasons)) {

-        warning("MuTect stats file lacks failure_reasons column.",

+        flog.warn("MuTect stats file lacks failure_reasons column.",

             " Keeping all variants listed in stats file.")

-        return(filterVcfBasic(vcf, tumor.id.in.vcf, verbose=verbose, ...))

+        return(filterVcfBasic(vcf, tumor.id.in.vcf, ...))

     }

     n <- nrow(vcf)

@@ -279,9 +279,9 @@ verbose=TRUE, ... ){

     ids <- sort(unique(unlist(sapply(ignore, grep, stats$failure_reasons))))

     vcf <- vcf[-ids]

-    if (verbose) message("Removing ", n-nrow(vcf), 

-        " MuTect calls due to blacklisted failure reasons.")

-    filterVcfBasic(vcf, tumor.id.in.vcf, verbose=verbose, ...)

+    flog.info("Removing %i MuTect calls due to blacklisted failure reasons.", 

+        n-nrow(vcf))

+    filterVcfBasic(vcf, tumor.id.in.vcf, ...)

 }

@@ -305,7 +305,6 @@ verbose=TRUE, ... ){

 #' value is for the case that variant is in both dbSNP and COSMIC > 2.

 #' @param tumor.id.in.vcf Id of tumor in case multiple samples are stored in

 #' VCF.

-#' @param verbose Verbose output.

 #' @return A \code{numeric(nrow(vcf))} vector with the prior probability of

 #' somatic status for each variant in the \code{CollapsedVCF}.

 #' @author Markus Riester

@@ -320,7 +319,7 @@ verbose=TRUE, ... ){

 #' @export setPriorVcf

 setPriorVcf <- function(vcf,

 prior.somatic=c(0.5, 0.0005, 0.999, 0.0001, 0.995, 0.01), 

-tumor.id.in.vcf=NULL, verbose=TRUE) {

+tumor.id.in.vcf=NULL) {

     if (is.null(tumor.id.in.vcf)) {

         tumor.id.in.vcf <- .getTumorIdInVcf(vcf)

     }

@@ -329,24 +328,24 @@ tumor.id.in.vcf=NULL, verbose=TRUE) {

          prior.somatic <- ifelse(info(vcf)$SOMATIC,

             prior.somatic[3],prior.somatic[4])

-         if (verbose) message("Found SOMATIC annotation in VCF. ",

-            "Setting somatic prior probabilities for somatic variants to ", 

-            tmp[3]," or to ", tmp[4], " otherwise.")

+         flog.info("Found SOMATIC annotation in VCF.")

+         flog.info("Setting somatic prior probabilities for somatic variants to %f or to %f otherwise.", 

+            tmp[3], tmp[4])

     } else {

          tmp <- prior.somatic

          prior.somatic <- ifelse(info(vcf)$DB,

             prior.somatic[2], prior.somatic[1])

          if (!is.null(info(vcf)$Cosmic.CNT)) {

-             if (verbose) message("Found COSMIC annotation in VCF. ",

-                "Setting somatic prior probabilities for hits to\n", tmp[5],

-                " or to ", tmp[6], " if in both COSMIC and dbSNP.")

+             flog.info("Found COSMIC annotation in VCF.")

+             flog.info("Setting somatic prior probabilities for hits to %f or to %f if in both COSMIC and dbSNP.", 

+                tmp[5], tmp[6])

              prior.somatic[which(info(vcf)$Cosmic.CNT>2)] <- tmp[5]

              prior.somatic[which(info(vcf)$Cosmic.CNT>2 & 

                 info(vcf)$DB)] <- tmp[6]

          } else {

-             if (verbose) message("Setting somatic prior probabilities ",

-                "for dbSNP hits to ", tmp[2]," or to ", tmp[1], " otherwise.")

+             flog.info("Setting somatic prior probabilities for dbSNP hits to %f or to %f otherwise.", 

+                tmp[2], tmp[1])

          }      

     }     

     prior.somatic

@@ -452,7 +451,7 @@ function(vcf, tumor.id.in.vcf, allowed=0.05) {

     }    

     return(TRUE)

 }    

-.readAndCheckVcf <- function(vcf.file, genome, verbose) {

+.readAndCheckVcf <- function(vcf.file, genome) {

     if (class(vcf.file) == "character") {    

         vcf <- readVcf(vcf.file, genome)

     } else if (class(vcf.file) != "CollapsedVCF") {

@@ -465,7 +464,7 @@ function(vcf, tumor.id.in.vcf, allowed=0.05) {

     if (sum(triAllelic)) {

         n <- nrow(vcf)

         vcf <- vcf[which(!triAllelic)]

-        if (verbose) message("Removing ",n-nrow(vcf), " triallelic sites.")

+        flog.info("Removing %i triallelic sites.",n-nrow(vcf))

     }    

     if (is.null(info(vcf)$DB)) {

         # try to add an DB field based on rownames

@@ -525,9 +524,9 @@ function(vcf, tumor.id.in.vcf, allowed=0.05) {

     vcf

 }

-.addCosmicCNT <- function(vcf, cosmic.vcf.file, verbose=TRUE) {

+.addCosmicCNT <- function(vcf, cosmic.vcf.file) {

     if (!is.null(info(vcf)$Cosmic.CNT)) {

-        if (verbose) message("VCF already COSMIC annotated. Skipping.")

+        flog.info("VCF already COSMIC annotated. Skipping.")

         return(vcf)        

     }

     cosmicSeqStyle <- seqlevelsStyle(headerTabix(

@@ -537,7 +536,7 @@ function(vcf, tumor.id.in.vcf, allowed=0.05) {

     if (!length(intersect(seqlevelsStyle(vcf), cosmicSeqStyle))) {

         seqlevelsStyle(vcfRenamedSL) <- cosmicSeqStyle[1]

     }

-    if (verbose) message("Reading COSMIC VCF...")

+    flog.info("Reading COSMIC VCF...")

     # look-up the variants in COSMIC

     cosmic.vcf <- readVcf(cosmic.vcf.file, genome=genome(vcf)[1],  

         ScanVcfParam(which = rowRanges(vcfRenamedSL),

@@ -569,20 +568,18 @@ function(vcf, tumor.id.in.vcf, allowed=0.05) {

     vcf

 }

-.annotateVcfTarget <- function(vcf, target.granges, interval.padding, verbose) {

+.annotateVcfTarget <- function(vcf, target.granges, interval.padding) {

     target.granges.padding <- target.granges

     start(target.granges.padding) <- start(target.granges.padding)-interval.padding

     end(target.granges.padding) <- end(target.granges.padding)+interval.padding

     .calcTargetedGenome <- function(granges) {

         tmp <- reduce(granges)

-        round(sum(width(tmp))/(1000^2),digits=2)

-    }    

-    if (verbose) {

-        message("Total size of targeted genomic region: ", 

-            .calcTargetedGenome(target.granges), "Mb (",

-            .calcTargetedGenome(target.granges.padding),

-            "Mb with ", interval.padding, "bp padding)")

+        sum(width(tmp))/(1000^2)