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

 .filterVcfByBQ <- function(vcf, tumor.id.in.vcf, min.base.quality) {

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

     idx <- .getBQFromVcf(vcf, tumor.id.in.vcf, base.quality.offset = 0) < min.base.quality

+    if (sum(idx, na.rm = TRUE) / length(idx) > 0.5) {

+        flog.warn("Many variants removed by min.base.quality (%i) filter. You might want to lower the cutoff.", min.base.quality)

+    }

     vcf <- .removeVariants(vcf, idx, "BQ", na.rm = FALSE)

     flog.info("Removing %i low quality variants with non-offset BQ < %i.",

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

@@ -42,6 +42,9 @@ option_list <- list(

     make_option(c("--base-quality-offset"), action = "store", type = "integer",

         default = formals(PureCN::filterVcfBasic)$base.quality.offset,

         help = "VCF Filter: Subtracts the specified value from the base quality score [default %default]"),

+    make_option(c("--min-base-quality"), action = "store", type = "integer",

+        default = formals(PureCN::filterVcfBasic)$min.base.quality,

+        help = "VCF Filter: Minimum base quality for variants. Set to 0 to turn off filter [default %default]"),

     make_option(c("--min-supporting-reads"), action = "store", type = "integer",

         default = formals(PureCN::filterVcfBasic)$min.supporting.reads,

         help = "VCF Filter: Instead of calculating the min. number of supporting reads, use specified one [default %default]"),

@@ -177,14 +180,14 @@ alias_list <- list(

     "maxhomozygousloss" = "max-homozygous-loss",

     "speedupheuristics" = "speedup-heuristics",

     "outvcf" = "out-vcf"

-)    

+)

 replace_alias <- function(x, deprecated = TRUE) {

     idx <- match(x, paste0("--", names(alias_list)))

     if (any(!is.na(idx))) {

         replaced <- paste0("--", alias_list[na.omit(idx)])

         x[!is.na(idx)] <- replaced

         if (deprecated) {

-            flog.warn("Deprecated arguments, use %s instead.", paste(replaced, collapse=" "))

+            flog.warn("Deprecated arguments, use %s instead.", paste(replaced, collapse = " "))

         }

     }

     return(x)

@@ -253,9 +256,9 @@ if (Sys.getenv("PURECN_DEBUG") != "") {

 }

 .checkFileList <- function(file) {

-    files <- read.delim(file, as.is = TRUE, header = FALSE)[,1]

+    files <- read.delim(file, as.is = TRUE, header = FALSE)[, 1]

     numExists <- sum(file.exists(files), na.rm = TRUE)

-    if (numExists < length(files)) { 

+    if (numExists < length(files)) {

         stop("File not exists in file ", file)

     }

     files

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

                 af.range = af.range, stats.file = opt$stats_file,

                 min.supporting.reads = opt$min_supporting_reads,

                 base.quality.offset = opt$base_quality_offset,

+                min.base.quality = opt$min_base_quality,

                 ignore = mutect.ignore,

                 interval.padding = opt$interval_padding),

             args.filterIntervals = list(

@@ -133,3 +133,10 @@ test_that("issue 249", {

     vcf.bq <- PureCN:::.readAndCheckVcf(vcf.bq, "hg19")

     expect_equal(NA, geno(vcf.bq)$BQ[[2, 1]])

 })

+

+test_that("issue 320", {

+    # set random BQ to NA

+    expect_output(x <- filterVcfBasic(vcf, use.somatic.status = FALSE, min.base.quality = 34),

+        "Many variants removed by min.base.quality")

+    expect_equal(34, min(sapply(geno(x$vcf)$BQ[,1], function(x) x[1])))

+})

@@ -207,7 +207,10 @@ Important recommendations:

   to presence in germline databases. _Mutect 1.1.7_ automatically calls SNPs, 

   but _Mutect 2_ does not. Make sure to run _Mutect 2_ with

   `--genotype-germline-sites true --genotype-pon-sites true`. You will not get

-  usuable output without those flags.

+  usuable output without those flags. Since _Mutect 2_ from _GATK 4.2.0+_,

+  average base quality scores can be very low and variants will be too

+  aggressively removed by _PureCN_.  You will need to set `--min-base-quality

+  20` in _PureCN.R_ to keep them.

 - Run the variant caller with a 50-75 base pair interval padding to increase

   the number of heterozygous SNPs (for example `--interval_padding` and

@@ -524,7 +527,9 @@ Important recommendations:

       with 50-100bp interval padding or no interval file at all. Also check

       that the interval file was generated using the baits coordinates, not the

       targets (the baits BED file should have a more even size distribution,

-      e.g. 120bp and multiples of it).

+      e.g. 120bp and multiples of it). If many variants are removed by the

+      default 25 base quality feature, you might be using _Mutect 2_ and need

+      to re-run _PureCN.R_ with `--min-base-quality 20`.

     - "Initial testing for significant sample cross-contamination" in the log

       file should not have many false positives, i.e. should be "unlikely" for