From 5169999d03b10c9715039bd842b42bea73cef526 Mon Sep 17 00:00:00 2001 From: Markus Riester Date: Sun, 31 Mar 2024 16:38:37 -0400 Subject: [PATCH 1/3] Attempt in a resolution for #40. --- DESCRIPTION | 2 +- NAMESPACE | 1 + NEWS | 6 ++++ R/adjustLogRatio.R | 40 +++++++++++++++++++++++ man/adjustLogRatio.Rd | 47 ++++++++++++++++++++++++++++ man/filterVcfMuTect2.Rd | 6 ++-- tests/testthat/test_adjustLogRatio.R | 15 +++++++++ 7 files changed, 113 insertions(+), 4 deletions(-) create mode 100644 R/adjustLogRatio.R create mode 100644 man/adjustLogRatio.Rd create mode 100644 tests/testthat/test_adjustLogRatio.R diff --git a/DESCRIPTION b/DESCRIPTION index 3d312009..6b8dc285 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -69,4 +69,4 @@ biocViews: CopyNumberVariation, Software, Sequencing, VariantAnnotation, VariantDetection, Coverage, ImmunoOncology NeedsCompilation: no ByteCompile: yes -RoxygenNote: 7.2.3.9000 +RoxygenNote: 7.3.1 diff --git a/NAMESPACE b/NAMESPACE index ca9ea465..cc7f1f3e 100755 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,5 +1,6 @@ # Generated by roxygen2: do not edit by hand +export(adjustLogRatio) export(annotateTargets) export(bootstrapResults) export(calculateBamCoverageByInterval) diff --git a/NEWS b/NEWS index 6ecd9d48..9405bd5d 100755 --- a/NEWS +++ b/NEWS @@ -1,6 +1,12 @@ Changes in version 2.10.0 ------------------------- +NEW FEATURES + o adjustLogRatio function for adjusting a tumor vs normal coverage + ratio for purity and ploidy. Useful for downstream tools that + expect ratios instead of absolute copy numbers such as GISTIC. + Thanks @tedtoal (#40). + SIGNIFICANT USER-VISIBLE CHANGES o Provide interval-level likelihood scores in runAbsoluteCN return diff --git a/R/adjustLogRatio.R b/R/adjustLogRatio.R new file mode 100644 index 00000000..c07c0383 --- /dev/null +++ b/R/adjustLogRatio.R @@ -0,0 +1,40 @@ +#' Adjust tumor vs. normal coverage log ratio for tumor purity and ploidy +#' +#' This function can be used to adjust the log ratio for tumor purity and +#' ploidy for downstream tools that expect a log2 ratio (for example GISTIC). +#' +#' +#' @param ratio Vector of log2 tumor vs normal coverage ratios. +#' @param purity Purity of sample. +#' @param ploidy Ploidy of sample. +#' @param is.log2 \code{log.ratio} is \code{log2} transformed. +#' @param min.ratio Minimum (non-log2-transformed) ratio. Set to approx -8 +#' \code{log2} adjusted. +#' @return \code{numeric(length(log.ratio))}, \code{log.ratio} adjusted +#' for \code{purity} and \code{ploidy} +#' @author Markus Riester +#' @references +# * Zack et al. (2012), Pan-cancer patterns of somatic copy number alteration +#' Nature Biotechnology. +#' * Toal (2018), https://github.com/lima1/PureCN/issues/40 +#' +#' @examples +#' +#' normal.coverage.file <- system.file("extdata", "example_normal.txt.gz", +#' package = "PureCN") +#' tumor.coverage.file <- system.file("extdata", "example_tumor.txt.gz", +#' package = "PureCN") +#' normal <- readCoverageFile(normal.coverage.file) +#' tumor <- readCoverageFile(tumor.coverage.file) +#' log.ratio <- calculateLogRatio(normal, tumor) +#' log.ratio.adjusted <- adjustLogRatio(log.ratio, 0.65, 1.73) +#' +#' @export adjustLogRatio +adjustLogRatio <- function(ratio, purity, ploidy, is.log2 = TRUE, min.ratio = 0.004) { + if (is.log2) ratio <- 2^ratio + adjusted <- (purity * ploidy * ratio + 2 * (1 - purity) * ratio - 2 * (1 - purity)) / (purity * ploidy) + adjusted <- pmax(min.ratio, adjusted) + if (is.log2) adjusted <- log2(adjusted) + return(adjusted) +} + diff --git a/man/adjustLogRatio.Rd b/man/adjustLogRatio.Rd new file mode 100644 index 00000000..0fb25abe --- /dev/null +++ b/man/adjustLogRatio.Rd @@ -0,0 +1,47 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/adjustLogRatio.R +\name{adjustLogRatio} +\alias{adjustLogRatio} +\title{Adjust tumor vs. normal coverage log ratio for tumor purity and ploidy} +\usage{ +adjustLogRatio(ratio, purity, ploidy, is.log2 = TRUE, min.ratio = 0.004) +} +\arguments{ +\item{ratio}{Vector of log2 tumor vs normal coverage ratios.} + +\item{purity}{Purity of sample.} + +\item{ploidy}{Ploidy of sample.} + +\item{is.log2}{\code{log.ratio} is \code{log2} transformed.} + +\item{min.ratio}{Minimum (non-log2-transformed) ratio. Set to approx -8 +\code{log2} adjusted.} +} +\value{ +\code{numeric(length(log.ratio))}, \code{log.ratio} adjusted +for \code{purity} and \code{ploidy} +} +\description{ +This function can be used to adjust the log ratio for tumor purity and +ploidy for downstream tools that expect a log2 ratio (for example GISTIC). +} +\examples{ + +normal.coverage.file <- system.file("extdata", "example_normal.txt.gz", + package = "PureCN") +tumor.coverage.file <- system.file("extdata", "example_tumor.txt.gz", + package = "PureCN") +normal <- readCoverageFile(normal.coverage.file) +tumor <- readCoverageFile(tumor.coverage.file) +log.ratio <- calculateLogRatio(normal, tumor) +log.ratio.adjusted <- adjustLogRatio(log.ratio, 0.65, 1.73) + +} +\references{ +Nature Biotechnology. + * Toal (2018), https://github.com/lima1/PureCN/issues/40 +} +\author{ +Markus Riester +} diff --git a/man/filterVcfMuTect2.Rd b/man/filterVcfMuTect2.Rd index ab88998e..cdc54ea1 100644 --- a/man/filterVcfMuTect2.Rd +++ b/man/filterVcfMuTect2.Rd @@ -7,9 +7,9 @@ filterVcfMuTect2( vcf, tumor.id.in.vcf = NULL, - ignore = c("clustered_events", "t_lod", "str_contraction", "read_position", - "position", "fragment_length", "multiallelic", "clipping", "strand_artifact", - "strand_bias", "slippage", "weak_evidence", "orientation", "haplotype"), + ignore = c("clustered_events", "t_lod", "str_contraction", "read_position", "position", + "fragment_length", "multiallelic", "clipping", "strand_artifact", "strand_bias", + "slippage", "weak_evidence", "orientation", "haplotype"), ... ) } diff --git a/tests/testthat/test_adjustLogRatio.R b/tests/testthat/test_adjustLogRatio.R new file mode 100644 index 00000000..c8e6062a --- /dev/null +++ b/tests/testthat/test_adjustLogRatio.R @@ -0,0 +1,15 @@ +context("adjustLogRatio") + +test_that("Function returns expected values for example coverage", { + data(purecn.example.output) + log.ratio <- purecn.example.output$results[[1]]$seg$seg.mean + purity <- purecn.example.output$results[[1]]$purity + ploidy <- purecn.example.output$results[[1]]$ploidy + log.ratio.adjusted <- adjustLogRatio(log.ratio, purity, ploidy) + total.ploidy <- 1.73 + p <- 1 + opt.C <- (2^(log.ratio.adjusted + log.ratio.offset) * total.ploidy)/p - ((2 * (1 - p))/p) + expect_lt(abs(min(log.ratio.adjusted, na.rm=TRUE) - log2(0.004)), 0.001) + expect_lt(median(abs(opt.C - purecn.example.output$results[[1]]$seg$C)), 0.1) +}) + From daf175bf1f3df9f06274d98a9e8c9e7b627826be Mon Sep 17 00:00:00 2001 From: Markus Riester Date: Sun, 31 Mar 2024 17:02:36 -0400 Subject: [PATCH 2/3] Fix broken test. --- tests/testthat/test_adjustLogRatio.R | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/testthat/test_adjustLogRatio.R b/tests/testthat/test_adjustLogRatio.R index c8e6062a..8c76a0cc 100644 --- a/tests/testthat/test_adjustLogRatio.R +++ b/tests/testthat/test_adjustLogRatio.R @@ -8,6 +8,7 @@ test_that("Function returns expected values for example coverage", { log.ratio.adjusted <- adjustLogRatio(log.ratio, purity, ploidy) total.ploidy <- 1.73 p <- 1 + log.ratio.offset <- 0 opt.C <- (2^(log.ratio.adjusted + log.ratio.offset) * total.ploidy)/p - ((2 * (1 - p))/p) expect_lt(abs(min(log.ratio.adjusted, na.rm=TRUE) - log2(0.004)), 0.001) expect_lt(median(abs(opt.C - purecn.example.output$results[[1]]$seg$C)), 0.1) From 5ec2226363f7e29b5886a87a5a32deb2a51fba88 Mon Sep 17 00:00:00 2001 From: Markus Riester Date: Mon, 1 Apr 2024 11:39:12 -0400 Subject: [PATCH 3/3] Changed default of min.ratio to cleaner -8. --- R/adjustLogRatio.R | 2 +- man/adjustLogRatio.Rd | 2 +- tests/testthat/test_adjustLogRatio.R | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/R/adjustLogRatio.R b/R/adjustLogRatio.R index c07c0383..2fe0770b 100644 --- a/R/adjustLogRatio.R +++ b/R/adjustLogRatio.R @@ -30,7 +30,7 @@ #' log.ratio.adjusted <- adjustLogRatio(log.ratio, 0.65, 1.73) #' #' @export adjustLogRatio -adjustLogRatio <- function(ratio, purity, ploidy, is.log2 = TRUE, min.ratio = 0.004) { +adjustLogRatio <- function(ratio, purity, ploidy, is.log2 = TRUE, min.ratio = 2^-8) { if (is.log2) ratio <- 2^ratio adjusted <- (purity * ploidy * ratio + 2 * (1 - purity) * ratio - 2 * (1 - purity)) / (purity * ploidy) adjusted <- pmax(min.ratio, adjusted) diff --git a/man/adjustLogRatio.Rd b/man/adjustLogRatio.Rd index 0fb25abe..fe70bc2b 100644 --- a/man/adjustLogRatio.Rd +++ b/man/adjustLogRatio.Rd @@ -4,7 +4,7 @@ \alias{adjustLogRatio} \title{Adjust tumor vs. normal coverage log ratio for tumor purity and ploidy} \usage{ -adjustLogRatio(ratio, purity, ploidy, is.log2 = TRUE, min.ratio = 0.004) +adjustLogRatio(ratio, purity, ploidy, is.log2 = TRUE, min.ratio = 2^-8) } \arguments{ \item{ratio}{Vector of log2 tumor vs normal coverage ratios.} diff --git a/tests/testthat/test_adjustLogRatio.R b/tests/testthat/test_adjustLogRatio.R index 8c76a0cc..830a257e 100644 --- a/tests/testthat/test_adjustLogRatio.R +++ b/tests/testthat/test_adjustLogRatio.R @@ -10,7 +10,7 @@ test_that("Function returns expected values for example coverage", { p <- 1 log.ratio.offset <- 0 opt.C <- (2^(log.ratio.adjusted + log.ratio.offset) * total.ploidy)/p - ((2 * (1 - p))/p) - expect_lt(abs(min(log.ratio.adjusted, na.rm=TRUE) - log2(0.004)), 0.001) + expect_lt(abs(min(log.ratio.adjusted, na.rm=TRUE) + 8), 0.001) expect_lt(median(abs(opt.C - purecn.example.output$results[[1]]$seg$C)), 0.1) })