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

 Type: Package

 Title: Copy number calling and SNV classification using

     targeted short read sequencing

-Version: 1.23.18

-Date: 2021-08-20

+Version: 1.23.19

+Date: 2021-08-21

 Authors@R: c(person("Markus", "Riester",

                     role = c("aut", "cre"),

                     email = "[email protected]",

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

     VGAM,

     tools,

     methods,

+    mclust,

     rhdf5,

     Matrix

 Suggests:

@@ -139,6 +139,9 @@ importFrom(graphics,strwidth)

 importFrom(graphics,symbols)

 importFrom(graphics,text)

 importFrom(gridExtra,grid.arrange)

+importFrom(mclust,Mclust)

+importFrom(mclust,emControl)

+importFrom(mclust,mclustBIC)

 importFrom(methods,is)

 importFrom(rhdf5,H5Fopen)

 importFrom(rtracklayer,import)

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

     o Report median absolute pairwise difference (MAPD) of tumor vs normal log2

       ratios in runAbsoluteCN

+    o Improved mapping bias estimates: variants with insufficient information

+      for position-specific fits (default 3-6 heterozygous variants) 

+      are clustered and assigned to the most similar fit  

 SIGNIFICANT USER-VISIBLE CHANGES

@@ -25,6 +28,7 @@ SIGNIFICANT USER-VISIBLE CHANGES

       finding noisy samples or batches in normal databases  

     o Do not error out readCurationFile when CSV is missing and directory

       is not writable when re-generating it (#196)

+    o Add segmentation parameters as attributes to segmentation data.frame  

 BUGFIXES

@@ -10,9 +10,12 @@

 #' @param min.normals Minimum number of normals with heterozygous SNP for

 #' calculating position-specific mapping bias.

 #' @param min.normals.betafit Minimum number of normals with heterozygous SNP

-#' fitting a beta distribution

+#' fitting a beta binomial distribution

+#' @param min.normals.assign.betafit Minimum number of normals with 

+#' heterozygous SNPs to assign to a beta binomal fit cluster

 #' @param min.median.coverage.betafit Minimum median coverage of normals with

 #' heterozygous SNP for fitting a beta distribution

+#' @param num.betafit.clusters Maximum number of beta binomial fit clusters

 #' @param yieldSize See \code{TabixFile}

 #' @param genome See \code{readVcf}

 #' @return A \code{GRanges} object with mapping bias and number of normal

@@ -28,12 +31,15 @@

 #' @importFrom GenomicRanges GRangesList

 #' @importFrom VGAM vglm Coef betabinomial dbetabinom

 #' @importFrom data.table rbindlist

+#' @importFrom mclust Mclust mclustBIC emControl

 #' @export calculateMappingBiasVcf

 calculateMappingBiasVcf <- function(normal.panel.vcf.file,

                                     min.normals = 1,

                                     min.normals.betafit = 7,

+                                    min.normals.assign.betafit = 3,

                                     min.median.coverage.betafit = 5,

-                                    yieldSize = 5000, genome) {

+                                    num.betafit.clusters = 9,

+                                    yieldSize = 50000, genome) {

     tab <- TabixFile(normal.panel.vcf.file, yieldSize = yieldSize)

     open(tab)

     param <- ScanVcfParam(geno = c("AD"), fixed = "ALT", info = NA)

@@ -48,7 +54,10 @@ calculateMappingBiasVcf <- function(normal.panel.vcf.file,

         mappingBias <- .calculateMappingBias(nvcf = vcf_yield,

             min.normals = min.normals,

             min.normals.betafit = min.normals.betafit,

-            min.median.coverage.betafit = min.median.coverage.betafit)

+            min.normals.assign.betafit = min.normals.assign.betafit,

+            min.median.coverage.betafit = min.median.coverage.betafit,

+            num.betafit.clusters = num.betafit.clusters

+        )

         if (length(ret)) {

             ret <- append(ret, GRangesList(mappingBias))

         } else {

@@ -61,7 +70,9 @@ calculateMappingBiasVcf <- function(normal.panel.vcf.file,

     attr(bias, "normal.panel.vcf.file") <- normal.panel.vcf.file

     attr(bias, "min.normals") <- min.normals

     attr(bias, "min.normals.betafit") <- min.normals.betafit

+    attr(bias, "min.normals.assign.betafit") <- min.normals.assign.betafit

     attr(bias, "min.median.coverage.betafit") <- min.median.coverage.betafit

+    attr(bias, "num.betafit.clusters") <- num.betafit.clusters

     attr(bias, "genome") <- genome

     bias

 }

@@ -78,8 +89,11 @@ calculateMappingBiasVcf <- function(normal.panel.vcf.file,

 #' calculating position-specific mapping bias.

 #' @param min.normals.betafit Minimum number of normals with heterozygous SNP

 #' fitting a beta distribution

+#' @param min.normals.assign.betafit Minimum number of normals with 

+#' heterozygous SNPs to assign to a beta binomal fit cluster

 #' @param min.median.coverage.betafit Minimum median coverage of normals with

 #' heterozygous SNP for fitting a beta distribution

+#' @param num.betafit.clusters Maximum number of beta binomial fit clusters

 #' @param AF.info.field Field in the \code{workspace} that stores the allelic 

 #' fraction

 #' @return A \code{GRanges} object with mapping bias and number of normal

@@ -104,7 +118,9 @@ calculateMappingBiasVcf <- function(normal.panel.vcf.file,

 calculateMappingBiasGatk4 <- function(workspace, reference.genome,

                                     min.normals = 1,

                                     min.normals.betafit = 7,

+                                    min.normals.assign.betafit = 3,

                                     min.median.coverage.betafit = 5,

+                                    num.betafit.clusters = 9,

                                     AF.info.field = "AF") {

     if (!requireNamespace("genomicsdb", quietly = TRUE) || 

@@ -160,7 +176,9 @@ calculateMappingBiasGatk4 <- function(workspace, reference.genome,

         alt = m_alt, ref = m_ref, gr = gr,

         min.normals = min.normals,

         min.normals.betafit = min.normals.betafit,

+        min.normals.assign.betafit = min.normals.assign.betafit,

         min.median.coverage.betafit = min.median.coverage.betafit,

+        num.betafit.clusters = num.betafit.clusters,

         verbose = TRUE

     )

     attr(bias, "workspace") <- workspace

@@ -187,7 +205,9 @@ calculateMappingBiasGatk4 <- function(workspace, reference.genome,

 .calculateMappingBias <- function(nvcf, alt = NULL, ref = NULL, gr = NULL,

                                   min.normals, min.normals.betafit = 7,

+                                  min.normals.assign.betafit = 3,

                                   min.median.coverage.betafit = 5,

+                                  num.betafit.clusters = 9,

                                   verbose = FALSE) {

     if (!is.null(nvcf)) {

         if (ncol(nvcf) < 2) {

@@ -240,6 +260,7 @@ calculateMappingBiasGatk4 <- function(workspace, reference.genome,

     gr$pon.count <- ponCntHits

     gr$mu <- x[5, ]

     gr$rho <- x[6, ]

+    gr <- .clusterFa(gr, fa, alt, ref, min.normals.assign.betafit, num.betafit.clusters)

     gr <- gr[order(gr$pon.count, decreasing = TRUE)]

     gr <- sort(gr)

     gr$triallelic <- FALSE

@@ -290,3 +311,27 @@ calculateMappingBiasGatk4 <- function(workspace, reference.genome,

     # get the alt allelic fraction for all SNPs

     apply(x, 2, function(y) y[2] / sum(head(y, 2)))

 }

+

+.clusterFa <- function(gr, fa, alt, ref, min.normals.assign.betafit = 3, num.betafit.clusters = 9) {

+    idx <- !is.na(gr$mu)

+    if (sum(idx) < num.betafit.clusters) return(gr) 

+    flog.info("Clustering beta binomial fits...")    

+    fit <- Mclust(mcols(gr)[idx, c("mu", "rho")], G = seq_len(num.betafit.clusters))

+    x <- sapply(seq_len(nrow(fa)), function(i) {

+        idx <- !is.na(fa[i, ]) & fa[i, ] > 0.05 & fa[i, ] < 0.9

+        if (!sum(idx) >= min.normals.assign.betafit) return(NA)

+        ll <- apply(fit$parameters$mean, 2, function(m) {

+                sum(dbetabinom(x = alt[i, idx], prob = m[1],

+                   size = alt[i, idx] + ref[i, idx],

+                   rho = m[2], log = TRUE))

+        })

+        if (is.infinite(max(ll))) return(NA)

+        which.max(ll)

+    })

+    idxa <- is.na(gr$mu) & !is.na(x)

+    flog.info("Assigning (%i/%i) variants a clustered beta binomal fit.",

+        sum(idxa), length(gr))

+    gr$mu[idxa] <- fit$parameters$mean[1,x[idxa]]

+    gr$rho[idxa] <- fit$parameters$mean[2,x[idxa]]

+    return(gr)

+}    

@@ -146,6 +146,14 @@ segmentationPSCBS <- function(normal, tumor, log.ratio, seg, plot.cnv,

             min.width = 3,alphaTCN = alpha, ...)

         try(flog.debug("Kappa: %f", PSCBS::estimateKappa(segPSCBS)), silent = TRUE)

         seg <- .PSCBSoutput2DNAcopy(segPSCBS, sampleid)

+        if (flog.threshold() == "DEBUG") {

+            flog.debug("Re-running segmentation without undo.SD...") 

+            segusd <- PSCBS::segmentByNonPairedPSCBS(input, tauA = tauA,

+                flavor = flavor, undoTCN = 0, knownSegments = knownSegments,

+                min.width = 3,alphaTCN = alpha, ...)

+            segusd <- .PSCBSoutput2DNAcopy(segusd, sampleid)

+            attr(seg, "Debug.undo.SD") <- segusd

+        }

         if (undo.SD <= 0 || is.null(max.segments) || nrow(seg) < max.segments) break

         flog.info("Found %i segments, exceeding max.segments threshold of %i.",

             nrow(seg), max.segments)

@@ -162,6 +170,17 @@ segmentationPSCBS <- function(normal, tumor, log.ratio, seg, plot.cnv,

     }

     seg <- .addAverageWeights(seg, weight.flag.pvalue, tumor, chr.hash)

     seg <- .fixBreakpointsInBaits(tumor, log.ratio, seg, chr.hash)

+    attr(seg, "PSCBS.Args") <- list(

+        alpha = alpha,

+        boost.on.target.max.size = boost.on.target.max.size,

+        flavor = flavor, 

+        iterations = try.again + 1,

+        min.logr.sdev = min.logr.sdev,

+        prune.hclust.h = prune.hclust.h,

+        prune.hclust.method = prune.hclust.method,

+        tauA = tauA,

+        undo.SD = undo.SD

+    )

     seg

 }

@@ -9,7 +9,9 @@ calculateMappingBiasGatk4(

   reference.genome,

   min.normals = 1,

   min.normals.betafit = 7,

+  min.normals.assign.betafit = 3,

   min.median.coverage.betafit = 5,

+  num.betafit.clusters = 9,

   AF.info.field = "AF"

 )

 }

@@ -24,9 +26,14 @@ calculating position-specific mapping bias.}

 \item{min.normals.betafit}{Minimum number of normals with heterozygous SNP

 fitting a beta distribution}

+\item{min.normals.assign.betafit}{Minimum number of normals with 

+heterozygous SNPs to assign to a beta binomal fit cluster}

+

 \item{min.median.coverage.betafit}{Minimum median coverage of normals with

 heterozygous SNP for fitting a beta distribution}

+\item{num.betafit.clusters}{Maximum number of beta binomial fit clusters}

+

 \item{AF.info.field}{Field in the \code{workspace} that stores the allelic 

 fraction}

 }

@@ -8,8 +8,10 @@ calculateMappingBiasVcf(

   normal.panel.vcf.file,

   min.normals = 1,

   min.normals.betafit = 7,

+  min.normals.assign.betafit = 3,

   min.median.coverage.betafit = 5,

-  yieldSize = 5000,

+  num.betafit.clusters = 9,

+  yieldSize = 50000,

   genome

 )

 }

@@ -22,11 +24,16 @@ Needs to be compressed and indexed with bgzip and tabix, respectively.}

 calculating position-specific mapping bias.}

 \item{min.normals.betafit}{Minimum number of normals with heterozygous SNP

-fitting a beta distribution}

+fitting a beta binomial distribution}

+

+\item{min.normals.assign.betafit}{Minimum number of normals with 

+heterozygous SNPs to assign to a beta binomal fit cluster}

 \item{min.median.coverage.betafit}{Minimum median coverage of normals with

 heterozygous SNP for fitting a beta distribution}

+\item{num.betafit.clusters}{Maximum number of beta binomial fit clusters}

+

 \item{yieldSize}{See \code{TabixFile}}

 \item{genome}{See \code{readVcf}}