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

 Type: Package

 Title: Copy number calling and SNV classification using

     targeted short read sequencing

-Version: 1.21.13

+Version: 1.21.14

 Date: 2021-05-02

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

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

@@ -36,6 +36,7 @@ export(readLogRatioFile)

 export(readSegmentationFile)

 export(runAbsoluteCN)

 export(segmentationCBS)

+export(segmentationGATK4)

 export(segmentationHclust)

 export(segmentationPSCBS)

 export(setMappingBiasVcf)

@@ -163,6 +164,7 @@ importFrom(stats,runif)

 importFrom(stats,t.test)

 importFrom(stats,weighted.mean)

 importFrom(tools,file_ext)

+importFrom(utils,compareVersion)

 importFrom(utils,data)

 importFrom(utils,head)

 importFrom(utils,object.size)

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

       minimizes noise while keeping the resolution as high as possible

     o Added support for GATK4 CollectAllelicCounts output as alternative

       to Mutect  

+    o Added segmentationGATK4 to used GATK4's segmentation function 

+      ModelSegments 

 SIGNIFICANT USER-VISIBLE CHANGES

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

     gr

 }

-.writeLogRatioFileGATK4 <- function(x, file) {

-    gr <- x$input$log.ratio

+.writeLogRatioFileGATK4 <- function(x, id = 1, file) {

+    gr <- x$log.ratio

     if (is.null(gr$log.ratio)) {

         .stopRuntimeError("log.ratio NULL in .writeLogRatioFileGATK4")

     }

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

         LOG2_COPY_RATIO = gr$log.ratio

     )

     con <- file(file, open = "w")

-    .writeGATKHeader(x$input$vcf, id = 1, con, "log-ratio")

+    .writeGATKHeader(x$vcf, id, con, "log-ratio")

     write.table(output, con, row.names = FALSE, quote = FALSE, sep = "\t")

     close(con)

     invisible(output)

@@ -73,5 +73,6 @@ readLogRatioFile <- function(file, format, zero = NULL) {

         sl <- seqlengths(vcf)

         writeLines(paste("@SQ", paste0("SN:",names(sl)), paste0("LN:", sl), sep="\t"), con)

     }

-    writeLines(paste("@RG", "ID:PureCN", paste0("SM:", samples(header(vcf))[id]), sep="\t"), con)

+    sampleid <- .getSampleIdFromVcf(vcf, id)

+    writeLines(paste("@RG", "ID:PureCN", paste0("SM:", sampleid), sep="\t"), con)

 }    

@@ -35,13 +35,17 @@ readSegmentationFile <- function(seg.file, sampleid, model.homozygous = FALSE,

     return("DNAcopy")

 }

 .convertSegGATK4 <- function(seg, sampleid) {

+    lr_col_id <- "LOG2_COPY_RATIO_POSTERIOR_50"

+    if (is.null(seg[[lr_col_id]])) {

+        lr_col_id <- "MEAN_LOG2_COPY_RATIO"

+    }

     data.frame(

         ID = sampleid,

         chrom = seg$CONTIG,

         loc.start = seg$START,

         loc.end = seg$END,

         num.mark = seg$NUM_POINTS_COPY_RATIO,

-        seg.mean = seg$LOG2_COPY_RATIO_POSTERIOR_50

+        seg.mean = seg[[lr_col_id]]

     )

 }    

 .checkSeg <- function(seg, sampleid, model.homozygous, verbose=TRUE) {

new file mode 100644

@@ -0,0 +1,117 @@

+#' GATK4 ModelSegments segmentation function

+#' 

+#' A wrapper for GATK4s ModelSegmentation function, useful when normalization

+#' is performed with other tools than GATK4, for example PureCN.

+#' This function is called via the

+#' \code{fun.segmentation} argument of \code{\link{runAbsoluteCN}}.  The

+#' arguments are passed via \code{args.segmentation}.

+#' 

+#' 

+#' @param normal Coverage data for normal sample. Ignored in this function.

+#' @param tumor Coverage data for tumor sample.

+#' @param log.ratio Copy number log-ratios, one for each exon in coverage file.

+#' @param seg If segmentation was provided by the user, this data structure

+#' will contain this segmentation. Useful for minimal segmentation functions.

+#' Otherwise PureCN will re-segment the data. This segmentation function

+#' ignores this user provided segmentation.

+#' @param vcf Optional \code{CollapsedVCF} object with germline allelic ratios.

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

+#' VCF.

+#' @param normal.id.in.vcf Id of normal in in VCF. Currently not used.

+#' @param prune.hclust.h Ignored in this function.

+#' @param prune.hclust.method Ignored in this function.

+#' @param additional.cmd.args \code{character(1)}. By default,

+#' \code{ModelSegments} is called with default parameters. Provide additional 

+#' arguments here.

+#' @param chr.hash Not needed here since \code{ModelSegments} does not

+#' require numbered chromosome names.

+#' @param ... Currently unused arguments provided to other segmentation

+#' functions.

+#' @return \code{data.frame} containing the segmentation.

+#' @author Markus Riester

+#'

+#' @seealso \code{\link{runAbsoluteCN}}

+#' @examples

+#' 

+#' normal.coverage.file <- system.file("extdata", "example_normal_tiny.txt", 

+#'     package="PureCN")

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

+#'     package="PureCN")

+#' vcf.file <- system.file("extdata", "example.vcf.gz",

+#'     package="PureCN")

+#' 

+#' # The max.candidate.solutions, max.ploidy and test.purity parameters are set to

+#' # non-default values to speed-up this example.  This is not a good idea for real

+#' # samples.

+#'\dontrun{

+#'  ret <-runAbsoluteCN(normal.coverage.file=normal.coverage.file, 

+#'      tumor.coverage.file=tumor.coverage.file, vcf.file=vcf.file, 

+#'      sampleid="Sample1",  genome="hg19",

+#'      fun.segmentation = segmentationGATK4, max.ploidy=4,

+#'      test.purity=seq(0.3,0.7,by=0.05), max.candidate.solutions=1)

+#'}

+#' 

+#' @importFrom utils compareVersion

+#' @export segmentationGATK4

+segmentationGATK4 <- function(normal, tumor, log.ratio, seg, 

+    vcf = NULL, tumor.id.in.vcf = 1, normal.id.in.vcf = NULL,

+    prune.hclust.h = NULL, prune.hclust.method = NULL,

+    chr.hash = NULL, ...) {

+

+    min.version <- "4.1.7.0"

+    if (.checkGATK4Version(min.version) < 0) {

+        .stopUserError("Unable to find a recent (>= %s) gatk binary.",

+            min.version)

+    }    

+    output.vcf.file <- NULL    

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

+        output.vcf.file <- tempfile(fileext = ".tsv")

+        .writeAllelicCountsFileGatk(vcf, tumor.id.in.vcf, output.vcf.file)

+    } 

+    if (length(log.ratio) != length(tumor)) {

+        .stopRuntimeError("tumor and log.ratio do not align in segmentationGATK4.")

+    }

+    output.lr.file <- tempfile(fileext = ".tsv")

+    # following .writeLogRatioFileGATK4 needs the GRanges information, so

+    # simply update it (should be the same though!)

+    tumor$log.ratio <- log.ratio

+    .writeLogRatioFileGATK4(list(vcf=vcf, log.ratio = tumor), tumor.id.in.vcf, output.lr.file)

+    output.dir <- tempdir()

+    sampleid <- .getSampleIdFromVcf(vcf, tumor.id.in.vcf)

+    args <- paste("ModelSegments",

+              "--allelic-counts", output.vcf.file,

+              "--denoised-copy-ratios", output.lr.file,

+              "--output-prefix", "tumor",

+              "-O", output.dir)

+    output <- try(system2("gatk", args, stderr = TRUE, stdout = TRUE))

+    if (is(output, "try-error")) {

+        .stopUserError("GATK4 ModelSegments runtime error: ",

+            output, "\nArguments: ", args)

+    }

+    seg <- readSegmentationFile(file.path(output.dir, "tumor.modelFinal.seg"), sampleid)

+    idx.enough.markers <- seg$num.mark > 1

+    rownames(seg) <- NULL

+    file.remove(output.vcf.file)

+    file.remove(output.lr.file)

+    seg[idx.enough.markers,]

+}

+

+.getSampleIdFromVcf <- function(vcf, tumor.id.in.vcf) {

+    if (is(tumor.id.in.vcf, "character") &&

+        tumor.id.in.vcf %in% samples(header(vcf))) return(tumor.id.in.vcf)

+    if (is(tumor.id.in.vcf, "numeric") &&

+        length(samples(header(vcf))) <= tumor.id.in.vcf) {

+        return(samples(header(vcf))[tumor.id.in.vcf])

+    }

+    return(NA)    

+}    

+

+.checkGATK4Version <- function(min.version) {

+    output <- try(system2("gatk", "--version", stdout=TRUE), silent = TRUE)

+    if (is(output, "try-error")) {

+        flog.warn("Cannot find gatk binary in path.")

+        return(-1)

+    }

+    prefix <- "The Genome Analysis Toolkit (GATK) v"    

+    compareVersion(gsub(prefix, "", output[grep(prefix, output, fixed = TRUE)], fixed = TRUE), min.version)

+}    

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

 #' \code{args.segmentation}.

 #' 

 #' 

-#' @param normal Coverage data for normal sample.

+#' @param normal Coverage data for normal sample. Ignored in this function.

 #' @param tumor Coverage data for tumor sample.

 #' @param log.ratio Copy number log-ratios, one for each exon in coverage file.

 #' @param seg If segmentation was provided by the user, this data structure

@@ -58,7 +58,7 @@ option_list <- list(

         default = formals(PureCN::filterIntervals)$min.total.counts,

         help = "Interval Filter: Keep only intervals with at least that many counts in both tumor and (tanget) normal [default %default]"),

     make_option(c("--funsegmentation"), action = "store", type = "character", default = "CBS",

-        help = "Segmentation: Algorithm. CBS, PSCBS, Hclust, or none [default %default]"),

+        help = "Segmentation: Algorithm. CBS, PSCBS, GATK4, Hclust, or none [default %default]"),

     make_option(c("--alpha"), action = "store", type = "double",

         default = formals(PureCN::segmentationCBS)$alpha,

         help = "Segmentation: significance of breakpoints [default %default]"),

@@ -260,6 +260,8 @@ if (file.exists(file.rds) && !opt$force) {

         } else if (opt$funsegmentation == "Hclust") {

             fun.segmentation <- segmentationHclust

             if (!is.null(seg.file)) uses.recommended.fun <- TRUE

+        } else if (opt$funsegmentation == "GATK4") {

+            fun.segmentation <- segmentationGATK4

         } else if (opt$funsegmentation == "none") {

             fun.segmentation <- function(seg, ...) seg

         } else {

@@ -297,7 +299,7 @@ if (file.exists(file.rds) && !opt$force) {

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

     }    

     vcf <- opt$vcf

-    if (!is.null(vcf) && file_ext(vcf) == "tsv") {

+    if (!is.null(vcf) && tools::file_ext(vcf) == "tsv") {

         flog.info("*.tsv file provided for --vcf, assuming GATK4 CollectAllelicCounts format")

         vcf <- readAllelicCountsFile(vcf)

     }    

new file mode 100644

@@ -0,0 +1,90 @@

+% Generated by roxygen2: do not edit by hand

+% Please edit documentation in R/segmentationGATK4.R

+\name{segmentationGATK4}

+\alias{segmentationGATK4}

+\title{GATK4 ModelSegments segmentation function}

+\usage{

+segmentationGATK4(

+  normal,

+  tumor,

+  log.ratio,

+  seg,

+  vcf = NULL,

+  tumor.id.in.vcf = 1,

+  normal.id.in.vcf = NULL,

+  prune.hclust.h = NULL,

+  prune.hclust.method = NULL,

+  chr.hash = NULL,

+  ...

+)

+}

+\arguments{

+\item{normal}{Coverage data for normal sample. Ignored in this function.}

+

+\item{tumor}{Coverage data for tumor sample.}

+

+\item{log.ratio}{Copy number log-ratios, one for each exon in coverage file.}

+

+\item{seg}{If segmentation was provided by the user, this data structure

+will contain this segmentation. Useful for minimal segmentation functions.

+Otherwise PureCN will re-segment the data. This segmentation function

+ignores this user provided segmentation.}

+

+\item{vcf}{Optional \code{CollapsedVCF} object with germline allelic ratios.}

+

+\item{tumor.id.in.vcf}{Id of tumor in case multiple samples are stored in

+VCF.}

+

+\item{normal.id.in.vcf}{Id of normal in in VCF. Currently not used.}

+

+\item{prune.hclust.h}{Ignored in this function.}

+

+\item{prune.hclust.method}{Ignored in this function.}

+

+\item{chr.hash}{Not needed here since \code{ModelSegments} does not

+require numbered chromosome names.}

+

+\item{...}{Currently unused arguments provided to other segmentation

+functions.}

+

+\item{additional.cmd.args}{\code{character(1)}. By default,

+\code{ModelSegments} is called with default parameters. Provide additional 

+arguments here.}

+}

+\value{

+\code{data.frame} containing the segmentation.

+}

+\description{

+A wrapper for GATK4s ModelSegmentation function, useful when normalization

+is performed with other tools than GATK4, for example PureCN.

+This function is called via the

+\code{fun.segmentation} argument of \code{\link{runAbsoluteCN}}.  The

+arguments are passed via \code{args.segmentation}.

+}

+\examples{

+

+normal.coverage.file <- system.file("extdata", "example_normal_tiny.txt", 

+    package="PureCN")

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

+    package="PureCN")

+vcf.file <- system.file("extdata", "example.vcf.gz",

+    package="PureCN")

+

+# The max.candidate.solutions, max.ploidy and test.purity parameters are set to

+# non-default values to speed-up this example.  This is not a good idea for real

+# samples.

+\dontrun{

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

+     tumor.coverage.file=tumor.coverage.file, vcf.file=vcf.file, 

+     sampleid="Sample1",  genome="hg19",

+     fun.segmentation = segmentationGATK4, max.ploidy=4,

+     test.purity=seq(0.3,0.7,by=0.05), max.candidate.solutions=1)

+}

+

+}

+\seealso{

+\code{\link{runAbsoluteCN}}

+}

+\author{

+Markus Riester

+}

@@ -29,7 +29,7 @@ segmentationPSCBS(

 )

 }

 \arguments{

-\item{normal}{Coverage data for normal sample.}

+\item{normal}{Coverage data for normal sample. Ignored in this function.}

 \item{tumor}{Coverage data for tumor sample.}

@@ -16,17 +16,17 @@ test_that("Example data matches", {

 })

 test_that("parsing -> writing -> parsing works", {

-    x <- purecn.example.output

+    x <- purecn.example.output$input

     y <- x

-    y$input$log.ratio$log.ratio <- NULL

+    y$log.ratio$log.ratio <- NULL

     output.file <- tempfile(fileext = ".tsv")

     expect_error(

-        PureCN:::.writeLogRatioFileGATK4(y, output.file),

+        PureCN:::.writeLogRatioFileGATK4(y, 1, output.file),

         "log.ratio NULL"      

     )

-    PureCN:::.writeLogRatioFileGATK4(x, output.file)

+    PureCN:::.writeLogRatioFileGATK4(x, 1, output.file)

     z <- readLogRatioFile(output.file)

-    expect_equivalent(x$input$log.ratio$log.ratio, z$log.ratio)

+    expect_equivalent(x$log.ratio$log.ratio, z$log.ratio)

     file.remove(output.file)

 })

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

 context("segmentation")

+normal.coverage.file <- system.file("extdata", "example_normal_tiny.txt", 

+    package="PureCN")

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

+    package="PureCN")

+vcf.file <- system.file("extdata", "example.vcf.gz",

+    package="PureCN")

+

 test_that("Precomputed boudaries are correct", {

     data(purecn.DNAcopy.bdry)

     alpha <- formals(segmentationCBS)$alpha

@@ -10,3 +17,18 @@ test_that("Precomputed boudaries are correct", {

     sbdry <- getbdry(eta, nperm, max.ones)

     expect_equal(purecn.DNAcopy.bdry, sbdry)

 })

+

+

+test_that("GATK4 wrapper works for example data.", {

+   skip_if_not(PureCN:::.checkGATK4Version("4.1.7.0") >= 0, "gatk binary > 4.1.7.0 required")

+ 

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

+        tumor.coverage.file = tumor.coverage.file, vcf.file = vcf.file, 

+        sampleid = "Sample1",  genome = "hg19",

+        fun.segmentation = segmentationGATK4, max.ploidy = 4,

+        test.purity = seq(0.3, 0.7, by = 0.05),

+        max.candidate.solutions = 1, plot.cnv = FALSE)

+    

+    expect_equal(0.65, ret$results[[1]]$purity, tolerance = 0.02)

+    expect_equal(1.62, ret$results[[1]]$ploidy, tolerance = 0.2)

+})