PureCN best practices



Update from previous stable versions

r Biocpkg("PureCN") is backward compatible with input generated by versions 1.16 and 1.18. For versions 1.8 to 1.14, please re-run NormalDB.R (see also below):

$ Rscript $PURECN/NormalDB.R --outdir $OUT_REF \
    --coveragefiles example_normal.list \
    --genome hg19 --normal_panel $NORMAL_PANEL --assay agilent_v6

For upgrades from version 1.6, we highly recommend starting from scratch following this tutorial.


For the command line scripts described in this tutorial, we will need to install r Biocpkg("PureCN") with suggested dependencies:

if (!requireNamespace("BiocManager", quietly = TRUE))
BiocManager::install("PureCN", dependencies = TRUE)

Alternatively, manually install the packages required by the command line scripts:

BiocManager::install(c("PureCN", "optparse", 
    "TxDb.Hsapiens.UCSC.hg19.knownGene", "org.Hs.eg.db"))

(Replace hg19 with your genome version).

To use the alternative and in many cases recommended r CRANpkg("PSCBS") segmentation:

# default PSCBS without support of interval weights

# patched PSCBS with support of interval weights
BiocManager::install("lima1/PSCBS", ref="add_dnacopy_weighting")

To call mutational signatures, install the GitHub version of the r CRANpkg("deconstructSigs") package:


For the experimental support of importing variant calls from GATK4 GenomicsDB, follow the installations instructions from GenomicsDB-R.

Prepare environment and assay-specific reference files

system.file("extdata", package="PureCN")
$ export PURECN="/path/to/PureCN/extdata"
$ Rscript $PURECN/PureCN.R --help
Usage: /path/to/PureCN/inst/extdata/PureCN.R [options] ...
# specify path where PureCN should store reference files
$ export OUT_REF="reference_files"
$ Rscript $PURECN/IntervalFile.R --infile baits_hg19.bed \ 
    --fasta hg19.fa --outfile $OUT_REF/baits_hg19_intervals.txt \
    --offtarget --genome hg19 \
    --export $OUT_REF/baits_optimized_hg19.bed \
    --mappability wgEncodeCrgMapabilityAlign100mer.bigWig \
    --reptiming wgEncodeUwRepliSeqK562WaveSignalRep1.bigWig

Internally, this script uses r Biocpkg("rtracklayer") to parse the infile. Make sure that the file format matches the file extension. See the r Biocpkg("rtracklayer") documentation for problems loading the file. Check that the genome version of the baits file matches the reference. Do not include chrM baits in case the capture kit includes some.

The --offtarget flag will include off-target reads. Including them is recommended except for Amplicon data.

The --genome version is needed to annotate exons with gene symbols. Use hg19/hg38 for human genomes, not b37/b38. You might get a warning that an annotation package is missing. For hg19, install r Biocpkg("TxDb.Hsapiens.UCSC.hg19.knownGene") in R.

The --export argument is optional. If provided, this script will store the modified intervals as BED file for example (again every r Biocpkg("rtracklayer") format is supported). This is useful when the coverages are calculated with third-party tools like GATK.

The --mappability argument should provide a r Biocpkg("rtracklayer") parsable file with a mappability score in the first meta data column. If provided, off-target regions will be restricted to regions specified in this file. On-target regions with low mappability will be excluded. For hg19, download the file from the UCSC website. Choose the kmer size that best fits your average mapped read length. For hg38, download recommended 76-kmer or 100-kmer mappability files through the courtesy of the Waldron lab from:

See the FAQ section of the main vignette for instruction how to generate such a file for other references.

Similarly, the --reptiming argument takes a replication timing score in the same format. If provided, GC-normalized and log-transformed coverage is tested for a linear relationship with this score and normalized accordingly. This is optional and provides only a minor benefit for coverage normalization, but can identify samples with high proliferation. Requires --offtarget to be useful.

Create VCF files

r Biocpkg("PureCN") does not ship with a variant caller. Use a third-party tool to generate a VCF for each sample.

Important recommendations:

Run PureCN with internal segmentation

The following describes r Biocpkg("PureCN") runs with internal copy number normalization and segmentation.

What you will need:


For each sample, tumor and normal, calculate GC-normalized coverages:

# Calculate and GC-normalize coverage from a BAM file 
$ Rscript $PURECN/Coverage.R --outdir $OUT/$SAMPLEID \ 
    --bam ${SAMPLEID}.bam \
    --intervals $OUT_REF/baits_hg19_intervals.txt

# GC-normalize coverage from a GATK DepthOfCoverage file
Rscript $PURECN/Coverage.R --outdir $OUT/$SAMPLEID \
    --coverage ${SAMPLEID}.coverage.sample_interval_summary \ 
    --intervals $OUT_REF/baits_hg19_intervals.txt

Similar to GATK, this script also takes a text file containing a list of BAM or coverage file names (one per line). The file extension must be .list:

# Calculate and GC-normalize coverage from a list of BAM files 
$ Rscript $PURECN/Coverage.R --outdir $OUT \ 
    --bam normals.list \
    --intervals $OUT_REF/baits_hg19_intervals.txt \
    --cores 4

Important recommendations:


To build a normal database for coverage normalization, copy the paths to all GC-normalized normal coverage files in a single text file, line-by-line:

ls -a normal*loess.txt.gz | cat > example_normal.list

# From already GC-normalized files
$ Rscript $PURECN/NormalDB.R --outdir $OUT_REF \
    --coveragefiles example_normal.list \
    --genome hg19 --assay agilent_v6

# When normal panel VCF is available (highly recommended for
# unmatched samples)
$ Rscript $PURECN/NormalDB.R --outdir $OUT_REF \
    --coveragefiles example_normal.list \
    --normal_panel $NORMAL_PANEL \
    --genome hg19 \
    --assay agilent_v6

# For a Mutect2/GATK4 normal panel GenomicsDB (experimental)
$ Rscript $PURECN/NormalDB.R --outdir $OUT_REF \
    --coveragefiles example_normal.list \
    --normal_panel $GENOMICSDB-WORKSPACE-PATH/pon_db \
    --genome hg19 \
    --assay agilent_v6

Important recommendations:


Now that the assay-specific files are created and all coverages calculated, we run PureCN.R to normalize, segment and determine purity and ploidy:


# Without a matched normal (minimal test run)
$ Rscript $PURECN/PureCN.R --out $OUT/$SAMPLEID \
    --tumor $OUT/$SAMPLEID/${SAMPLEID}_coverage_loess.txt.gz \
    --sampleid $SAMPLEID \
    --vcf ${SAMPLEID}_mutect.vcf \
    --normaldb $OUT_REF/normalDB_hg19.rds \
    --intervals $OUT_REF/baits_hg19_intervals.txt \
    --genome hg19 

# Production pipeline run
$ Rscript $PURECN/PureCN.R --out $OUT/$SAMPLEID \
    --tumor $OUT/$SAMPLEID/${SAMPLEID}_coverage_loess.txt.gz \
    --sampleid $SAMPLEID \
    --vcf ${SAMPLEID}_mutect.vcf \
    --statsfile ${SAMPLEID}_mutect_stats.txt \
    --normaldb $OUT_REF/normalDB_hg19.rds \
    --mappingbiasfile $OUT_REF/mapping_bias_hg19.rds \
    --intervals $OUT_REF/baits_hg19_intervals.txt \
    --snpblacklist hg19_simpleRepeats.bed \
    --genome hg19 \
    --force --postoptimize --seed 123

# With a matched normal (test run; for production pipelines we recommend the
# unmatched workflow described above)
$ Rscript $PURECN/PureCN.R --out $OUT/$SAMPLEID \
    --tumor $OUT/$SAMPLEID/${SAMPLEID}_coverage_loess.txt.gz \
    --normal $OUT/$SAMPLEID/${SAMPLEID_NORMAL}_coverage_loess.txt.gz \
    --sampleid $SAMPLEID \
    --vcf ${SAMPLEID}_mutect.vcf \
    --normaldb $OUT_REF/normalDB_hg19.rds \
    --intervals $OUT_REF/baits_hg19_intervals.txt \
    --genome hg19

# Recreate output after manual curation of ${SAMPLEID}.csv
$ Rscript $PURECN/PureCN.R --rds $OUT/$SAMPLEID/${SAMPLEID}.rds

Important recommendations:

Run PureCN with third-party segmentation

What you will need:

General usage

If you already have a segmentation from third-party tools (for example CNVkit, GATK4, EXCAVATOR2). For a minimal test run:

Rscript $PURECN/PureCN.R --out $OUT/$SAMPLEID  \
    --sampleid $SAMPLEID \
    --segfile $OUT/$SAMPLEID/${SAMPLEID}.cnvkit.seg \
    --vcf ${SAMPLEID}_mutect.vcf \
    --intervals $OUT_REF/baits_hg19_intervals.txt \
    --genome hg19 

See the main vignette for more details and file formats.

Recommended CNVkit usage

For a production pipeline run we provide again more information about the assay and genome. Here an CNVkit example:

# Recommended: Provide a normal panel VCF to remove mapping biases, pre-compute
# position-specific bias for much faster runtimes with large panels
# This needs to be done only once for each assay
Rscript $PURECN/NormalDB.R --outdir $OUT_REF --normal_panel $NORMAL_PANEL \
    --assay agilent_v6 --genome hg19 --force

# Export the segmentation in DNAcopy format
cnvkit.py export seg $OUT/$SAMPLEID/${SAMPLEID}_cnvkit.cns --enumerate-chroms \
    -o $OUT/$SAMPLEID/${SAMPLEID}_cnvkit.seg

# Run PureCN by providing the *.cnr and *.seg files 
Rscript $PURECN/PureCN.R --out $OUT/$SAMPLEID  \
    --sampleid $SAMPLEID \
    --tumor $OUT/$SAMPLEID/${SAMPLEID}_cnvkit.cnr \
    --segfile $OUT/$SAMPLEID/${SAMPLEID}_cnvkit.seg \
    --mappingbiasfile $OUT_REF/mapping_bias_agilent_v6_hg19.rds \
    --vcf ${SAMPLEID}_mutect.vcf \
    --statsfile ${SAMPLEID}_mutect_stats.txt \
    --snpblacklist hg19_simpleRepeats.bed \
    --genome hg19 \
    --funsegmentation Hclust \
    --force --postoptimize --seed 123

Important recommendations:

Recommended GATK4 usage

# Recommended: Provide a normal panel GenomicsDB to remove mapping
# biases, pre-compute position-specific bias for much faster runtimes
# with large panels. This needs to be done only once for each assay. 
Rscript $PURECN/NormalDB.R --outdir $OUT_REF \
    --normal_panel $GENOMICSDB-WORKSPACE-PATH/pon_db \
    --assay agilent_v6 --genome hg19 --force

Rscript $PURECN/PureCN.R --out $OUT/$SAMPLEID  \
    --sampleid $SAMPLEID \
    --tumor $OUT/$SAMPLEID/${SAMPLEID}.hdf5 \
    --logratiofile $OUT/$SAMPLEID/${SAMPLEID}.denoisedCR.tsv \
    --segfile $OUT/$SAMPLEID/${SAMPLEID}.modelFinal.seg \
    --mappingbiasfile $OUT_REF/mapping_bias_agilent_v6_hg19.rds \
    --vcf ${SAMPLEID}_mutect2_filtered.vcf \
    --snpblacklist hg19_simpleRepeats.bed \
    --genome hg19 \
    --funsegmentation Hclust \
    --force --postoptimize --seed 123


Dx.R provides copy number and mutation metrics commonly used as biomarkers, most importantly tumor mutational burden (TMB), chromosomal instability (CIN) and mutational signatures.

# Provide a BED file with callable regions, for examples obtained by
# GATK CallableLoci. Useful to calculate mutations per megabase and
# to exclude low quality regions.
grep CALLABLE ${SAMPLEID}_callable_status.bed > \ 

# Only count mutations in callable regions, also subtract what was
# ignored in PureCN.R via --snpblacklist, like simple repeats, from the
# mutation per megabase calculation
# Also search for the COSMIC mutation signatures
# (http://cancer.sanger.ac.uk/cosmic/signatures)
Rscript $PureCN/Dx.R --out $OUT/$SAMPLEID/$SAMPLEID \
    --rds $OUT/SAMPLEID/${SAMPLEID}.rds \
    --callable ${SAMPLEID}_callable_status_filtered.bed \
    --exclude hg19_simpleRepeats.bed \

# Restrict mutation burden calculation to coding sequences
Rscript $PureCN/FilterCallableLoci.R --genome hg19 \
    --infile ${SAMPLEID}_callable_status_filtered.bed \
    --outfile ${SAMPLEID}_callable_status_filtered_cds.bed \
    --exclude '^HLA'

Rscript $PureCN/Dx.R --out $OUT/$SAMPLEID/${SAMPLEID}_cds \
    --rds $OUT/SAMPLEID/${SAMPLEID}.rds \
    --callable ${SAMPLEID}_callable_status_filtered_cds.bed \
    --exclude hg19_simpleRepeats.bed

Important recommendations:


Argument name | Corresponding PureCN argument | PureCN function -----------------------|-------------------------------|---------------- --fasta | reference.file | preprocessIntervals --infile | interval.file | preprocessIntervals --offtarget | off.target | preprocessIntervals --targetwidth | average.target.width | preprocessIntervals --mintargetwidth | min.target.width | preprocessIntervals --smalltargets | small.targets | preprocessIntervals --offtargetwidth | average.off.target.width | preprocessIntervals --offtargetseqlevels | off.target.seqlevels | preprocessIntervals --mappability | mappability | preprocessIntervals --minmappability | min.mappability | preprocessIntervals --reptiming | reptiming | preprocessIntervals --reptimingwidth | average.reptiming.width | preprocessIntervals
--genome | txdb, org | annotateTargets --outfile | | --export | | rtracklayer::export
--version -v | | --force -f | | --help -h | |

: (#tab:intervalfile) IntervalFile

Argument name | Corresponding PureCN argument | PureCN function -------------------|-------------------------------|---------------- --bam | bam.file | calculateBamCoverageByInterval --bai | index.file | calculateBamCoverageByInterval --coverage | coverage.file | correctCoverageBias --intervals | interval.file | correctCoverageBias --method | method | correctCoverageBias --keepduplicates | keep.duplicates | calculateBamCoverageByInterval --removemapq0 | mapqFilter | ScanBamParam --outdir | | --cores | | Number of CPUs to use when multiple BAMs are provided --parallel | | Use default r Biocpkg("BiocParallel") backend when multiple BAMs are provided --seed | | --version -v | | --force -f | | --help -h | |

: (#tab:coverage) Coverage

Argument name | Corresponding PureCN argument | PureCN function -----------------------|-------------------------------|---------------- --coveragefiles | normal.coverage.files | createNormalDatabase --normal_panel | normal.panel.vcf.file | calculateMappingBiasVcf --assay -a | Optional assay name | Used in output file names. --genome -g | Optional genome version | Used in output file names. --outdir -o | | --version -v | | --force -f | | --help -h | |

: (#tab:normaldb) NormalDB

Argument name | Corresponding PureCN argument | PureCN function -----------------------|-------------------------------|---------------- --sampleid -i | sampleid | runAbsoluteCN
--normal | normal.coverage.file | runAbsoluteCN
--tumor | tumor.coverage.file | runAbsoluteCN --vcf | vcf.file | runAbsoluteCN
--rds | file.rds | readCurationFile
--mappingbiasfile | mapping.bias.file | setMappingBiasVcf --normaldb | normalDB (serialized with saveRDS) | calculateTangentNormal, filterTargets --segfile | seg.file | runAbsoluteCN
--logratiofile | log.ratio | runAbsoluteCN
--additionaltumors | tumor.coverage.files | processMultipleSamples
--sex | sex | runAbsoluteCN
--genome | genome | runAbsoluteCN
--intervals | interval.file | runAbsoluteCN
--statsfile | stats.file | filterVcfMuTect
--minaf | af.range | filterVcfBasic --snpblacklist | snp.blacklist | filterVcfBasic
--error | error | runAbsoluteCN
--dbinfoflag | DB.info.flag | runAbsoluteCN
--popafinfofield | POPAF.info.field | runAbsoluteCN
--mincosmiccnt | min.cosmic.cnt | setPriorVcf
--funsegmentation | fun.segmentation | runAbsoluteCN
--alpha | alpha | segmentationCBS --undosd | undo.SD | segmentationCBS --maxsegments | max.segments | runAbsoluteCN --minpurity | test.purity | runAbsoluteCN --maxpurity | test.purity | runAbsoluteCN --minploidy | min.ploidy | runAbsoluteCN --maxploidy | max.ploidy | runAbsoluteCN --maxcopynumber | test.num.copy | runAbsoluteCN --postoptimize | post.optimize | runAbsoluteCN
--bootstrapn | n | bootstrapResults --modelhomozygous | model.homozygous | runAbsoluteCN
--model | model | runAbsoluteCN
--logratiocalibration | log.ratio.calibration | runAbsoluteCN
--maxnonclonal | max.non.clonal | runAbsoluteCN
--maxhomozygousloss | max.homozygous.loss | runAbsoluteCN
--outvcf | return.vcf | predictSomatic --out -o | |
--parallel | BPPARAM | runAbsoluteCN --cores | BPPARAM | runAbsoluteCN --seed | | --version -v | | --force -f | | --help -h | |

: (#tab:purecn) PureCN

Argument name | Corresponding PureCN argument | PureCN function ----------------|-------------------------------|---------------- --rds | file.rds | readCurationFile
--callable | callable | callMutationBurden --exclude | exclude | callMutationBurden --maxpriorsomatic | max.prior.somatic | callMutationBurden --signatures | | deconstructSigs::whichSignatures --signature_databases | | deconstructSigs::whichSignatures --out | | --version -v | | --force -f | | --help -h | |

: (#tab:dx) Dx

Try the PureCN package in your browser

Any scripts or data that you put into this service are public.

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