calculatePowerDetectSomatic: Power calculation for detecting somatic mutations

In PureCN: Copy number calling and SNV classification using targeted short read sequencing

Description

This function calculates the probability of correctly rejecting the null hypothesis that an alt allele is a sequencing error rather than a true (mono-)clonal mutation.

Usage

calculatePowerDetectSomatic(
  coverage,
  f = NULL,
  purity = NULL,
  ploidy = NULL,
  cell.fraction = 1,
  error = 0.001,
  fpr = 5e-07,
  verbose = TRUE
)

Arguments

coverage	Mean sequencing coverage.
f	Mean expected allelic fraction. If NULL, requires purity and ploidy and then calculates the expected fraction.
purity	Purity of sample. Only required when f is NULL.
ploidy	Ploidy of sample. Only required when f is NULL.
cell.fraction	Fraction of cells harboring mutation. Ignored if f is not NULL.
error	Estimated sequencing error rate.
fpr	Required false positive rate for mutation vs. sequencing error.
verbose	Verbose output.

Value

A list with elements

power	Power to detect somatic mutations.
k	Minimum number of supporting reads.
f	Expected allelic fraction.

Author(s)

Markus Riester

References

Carter et al. (2012), Absolute quantification of somatic DNA alterations in human cancer. Nature Biotechnology.

Examples

purity <- c(0.1,0.15,0.2,0.25,0.4,0.6,1)
coverage <- seq(5,35,1)
power <- lapply(purity, function(p) sapply(coverage, function(cv) 
    calculatePowerDetectSomatic(coverage=cv, purity=p, ploidy=2, 
    verbose=FALSE)$power))

# Figure S7b in Carter et al.
plot(coverage, power[[1]], col=1, xlab="Sequence coverage", 
    ylab="Detection power", ylim=c(0,1), type="l")

for (i in 2:length(power)) lines(coverage, power[[i]], col=i)
abline(h=0.8, lty=2, col="grey")     
legend("bottomright", legend=paste("Purity", purity), fill=seq_along(purity))

# Figure S7c in Carter et al.
coverage <- seq(5,350,1)
power <- lapply(purity, function(p) sapply(coverage, function(cv) 
    calculatePowerDetectSomatic(coverage=cv, purity=p, ploidy=2, 
        cell.fraction=0.2, verbose=FALSE)$power))
plot(coverage, power[[1]], col=1, xlab="Sequence coverage", 
    ylab="Detection power", ylim=c(0,1), type="l")

for (i in 2:length(power)) lines(coverage, power[[i]], col=i)
abline(h=0.8, lty=2, col="grey")     
legend("bottomright", legend=paste("Purity", purity), fill=seq_along(purity))

PureCN documentation built on Nov. 8, 2020, 5:37 p.m.