Bin size and resolution¶

If you need higher resolution, you can select a smaller average size for your target and antitarget bins.

Exons in the human genome have an average size of about 200bp. The target bin size default of 267 is chosen so that splitting larger exons will produce bins with a minimum size of 200. Since bins that contain fewer reads result in a noisier copy number signal, this approach ensures the “noisiness” of the bins produced by splitting larger exons will be no worse than average.

Setting the average size of target bins to 100bp, for example, will yield about twice as many target bins, which might result in higher-resolution segmentation. However, the number of reads counted in each bin will be reduced by about half, increasing the variance or “noise” in bin-level coverages. An excess of noisy bins can make visualization difficult, and since the noise may not be Gaussian, especially in the presence of many bins with zero reads, the CBS algorithm could produce less accurate segmentation results on low-coverage samples. In practice we see good results with an average of 200-300 reads per bin; we therefore recommend an overall on-target sequencing coverage depth of at least 200x to 300x with a read length of 100 to justify reducing the average target bin size to 100bp.

Adding gene names¶

In case the vendor BED file does not label each region with a corresponding gene name, the --annotate option can add or replace these labels. Gene annotation databases, e.g. RefSeq or Ensembl, are available in “flat” format from UCSC (e.g. refFlat.txt for hg19).

In other cases the region labels are a combination of human-readable gene names and database accession codes, separated by commas (e.g. “ref|BRAF,mRNA|AB529216,ens|ENST00000496384”). The --short-names option splits these accessions on commas, then chooses the single accession that covers in the maximum number of consecutive regions that share that accession, and applies it as the new label for those regions. (You may find it simpler to just apply the refFlat annotations.)

Off-target bin size¶

An appropriate off-target bin size can be computed as the product of the average target region size and the fold-enrichment of sequencing reads in targeted regions, such that roughly the same number of reads are mapped to on– and off-target bins on average — roughly proportional to the level of on-target enrichment.

The preliminary coverage information can be obtained with the script CalculateHsMetrics in the Picard suite (http://picard.sourceforge.net/), or from the console output of the CNVkit coverage command when run on the target regions.

BAM file preparation¶

For best results, use an aligner such as BWA-MEM, with the option to mark secondary mappings of reads, and flag PCR duplicates with a program such as SAMBLASTER, SAMBAMBA, or the MarkDuplicates script in Picard tools, so that CNVkit will skip these reads when calculating read depth.

You will probably want to index the finished BAM file using samtools or SAMBAMBA. But if you haven’t done this beforehand, CNVkit will automatically do it for you.

Paired or pooled normals¶

To analyze a cohort sequenced on a single platform, we recommend combining all normal samples into a pooled reference, even if matched tumor-normal pairs were sequenced – our benchmarking showed that a pooled reference performed slightly better than constructing a separate reference for each matched tumor-normal pair. Furthermore, even matched normals from a cohort sequenced together can exhibit distinctly different copy number biases (see Plagnol et al. 2012 and Backenroth et al. 2014); reusing a pooled reference across the cohort provides some consistency to help diagnose such issues.

Notes on sample selection:

• You can use cnvkit.py metrics *.cnr -s *.cns to see if any samples are especially noisy. See the metrics command.
• CNVkit will usually call larger CNAs reliably down to about 10x on-target coverage, but there will tend to be more spurious segments, and smaller-scale or subclonal CNAs can be hard to infer below that point. This is well below the minimum coverage thresholds typically used for SNV calling, especially for targeted sequencing of tumor samples that may have significant normal-cell contamination and subclonal tumor-cell populations. So, a normal sample that passes your other QC checks will probably be OK to use in building a CNVkit reference – assuming it was sequenced on the same platform as the other samples you’re calling.

If normal samples are not available, it will sometimes be acceptable to build the reference from a collection of tumor samples. You can use the scatter command on the raw .cnn coverage files to help choose samples with relatively minimal and non-recurrent CNVs for use in the reference.

With no control samples¶

Alternatively, you can create a “flat” reference of neutral copy number (i.e. log2 0.0) for each probe from the target and antitarget interval files. This still computes the GC content of each region if the reference genome is given.

cnvkit.py reference -o FlatReference.cnn -f ucsc.hg19.fa -t Tiled.bed -a Background.bed

Possible uses for a flat reference include:

1. Extract copy number information from one or a small number of tumor samples when no suitable reference or set of normal samples is available. The copy number calls will not be quite as accurate, but large-scale CNVs should still be visible.
2. Create a “dummy” reference to use as input to the batch command to process a set of normal samples. Then, create a “real” reference from the resulting *.targetcoverage.cnn and *.antitargetcoverage.cnn files, and re-run batch on a set of tumor samples using this updated reference.
3. Evaluate whether a given paired or pooled reference is suitable for an analysis by repeating the CNVkit analysis with a flat reference and comparing the CNAs found with both the original and flat reference for the same samples.

How it works¶

CNVkit uses robust methods to extract a usable signal from the reference samples.

Each input sample is first median-centered, then read-depth bias corrections (the same used in the fix command) are performed on each of the normal samples separately.

The samples’ median-centered, bias-corrected log2 read depth values are then combined to take the weighted average (Tukey’s biweight location) and spread (Tukey’s biweight midvariance) of the values at each on– and off-target genomic bin among all samples. (For background on these statistical methods see Lax (1985) and Randal (2008).) To adjust for the lower statistical reliability of a smaller number of samples for estimating parameters, a “pseudocount” equivalent to one sample of neutral copy number is included in the dataset when calculating these values.

These values are saved in the output “reference.cnn” file as the “log2” and “spread” columns, indicating the expected read depth and the reliability of this estimate.

If a FASTA file of the reference genome is given, for each genomic bin the fraction of GC (proportion of “G” and “C” characters among all “A”, “T”, “G” and “C” characters in the subsequence, ignoring “N” and any other ambiguous characters) and repeat-masked values (proportion of lowercased non-“N” characters in the sequence) are calculated and stored in the output “reference.cnn” file as columns “gc” and “rmask”. For efficiency, the samtools FASTA index file (.fai) is used to locate the binned sequence regions in the FASTA file. If the GC or RepeatMasker bias corrections are skipped using the --no-gc or --no-rmask options, then those columns are omitted from the output file; if both are skipped, then the genome FASTA file (if provided) is not examined at all.

The result is a reference copy-number profile that can then be used to correct other individual samples.

How it works¶

The “observed” on- and off-target read depths are each median-centered and bias-corrected, as when constructing the reference. The corresponding “expected” normalized log2 read-depth values from the reference are then subtracted for each set of bins.

Bias corrections use the GC and RepeatMasker information from the “gc” and “rmask” columns of the reference .cnn file; if those are missing (i.e. the reference was built without those corrections), fix will skip them too (with a warning). If you constructed the reference but then called fix with a different set of bias correction flags, the biases could be over- or under-corrected in the test sample – so use the options --no-gc, --no-rmask and --no-edge consistently or not at all.

CNVkit filters out bins failing certain predefined criteria: those where the reference log2 read depth is below a threshold (default -5), or the spread of read depths among all normal samples in the reference is above a threshold (default 1.0).

A weight is assigned to each remaining bin depending on:

1. The size of the bin;
2. The deviation of the bin’s log2 value in the reference from 0;
3. The “spread” of the bin in the reference.

(The latter two only apply if at least one normal/control sample was used to build the reference.)

Finally, the corrected on- and off-target bin-level copy ratios with associated weights are concatenated, sorted, and written to a .cnr file.

Transformations¶

With the --purity option, log2 ratios are rescaled using the same calculation as the rescale command. The observed log2 ratios in the input .cns file are treated as a mix of some fraction of tumor cells (specified by --purity), possibly with altered copy number, and a remainder of normal cells with neutral copy number (specified by --ploidy for autosomes; by default, diploid autosomes, haploid Y or X/Y depending on reference gender). This equation is rearranged to find the absolute copy number of the tumor cells alone, rounded to the nearest integer. The expected and observed ploidy of the sex chromosomes (X and Y) is different, so it’s important to specify -y/--male-reference if a male reference was used; the sample gender can be specified if known, otherwise it will be guessed from the average log2 ratio of chromosome X.

When a VCF file containing SNV calls for the same tumor sample (and optionally a matched normal) is given using the -v/--vcf option, the b-allele frequencies (BAFs) of the heterozygous, non-somatic SNVs falling within each segment are mirrored, averaged, and listed in the output .cns file as an additional “baf” column (using the same logic as export nexus-ogt). If --purity was specified, then the BAF values are also rescaled.

Calling methods¶

After the above adjustments, the “threshold” and “clonal” methods calculate the absolute integer copy number of each segment.

The “clonal” method converts the log2 values to absolute scale using the given --ploidy, then simply rounds the absolute copy number values to the nearest integer. This method is reasonable for germline samples, highly pure tumor samples (e.g. cell lines), or when the tumor fraction is accurately known and specified with --purity.

The “threshold” method applies fixed log2 ratio cutoff values for each integer copy number state. This approach can be an alternative to specifying and adjusting for the tumor cell fraction or purity directly. However, if --purity is specified, then the log2 values will still be rescaled before applying the copy number thresholds.

The default threshold values are reasonable for a tumor sample with purity of at least 40% or so. For germline samples, the -t values shown above may yield more accurate calls.

The thresholds map to integer copy numbers in order, starting from zero: log2 ratios up to the first threshold value are assigned a copy number 0, log2 ratios between the first and second threshold values get copy number 1, and so on.

For homogeneous samples of known ploidy, you can calculate cutoffs from scatch by log-transforming the integer copy number values of interest, plus .5 (for rounding), divided by the ploidy. For a diploid genome:

>>> import numpy as np
>>> copy_nums = np.arange(5)
>>> print(np.log2((copy_nums+.5) / 2)
[-2.         -0.4150375   0.32192809  0.80735492  1.169925  ]

Or, in R:

> log2( (0:4 + .5) / 2)
[1] -2.0000000 -0.4150375  0.3219281  0.8073549  1.1699250

Allele frequencies and counts¶

If a VCF file is given using the -v/--vcf option, then for each segment with a BAF value (i.e. where SNVs were available), allele-specific integer copy number values are inferred from the total copy number and BAF, and output in columns “cn1” and “cn2”. This calculation uses the same method as PSCBS: total copy number is multiplied by the BAF, and rounded to the nearest integer.

Allelic imbalance, including copy-number-neutral loss of heterozygosity (LOH), is then apparent when a segment’s “cn1” and “cn2” fields have different values.

bed¶

Segments can be exported to BED format to support a variety of other uses, such as viewing in a genome browser. The log2 ratio value of each segment is converted and rounded to an integer value, as required by the BED format. To get accurate copy number values, see the commands rescale and call.

# Estimate integer copy number of each segment
cnvkit.py call Sample.cns -y -o Sample.call.cns
# Show estimated integer copy number of all regions
cnvkit.py export bed Sample.call.cns --show all -y -o Sample.bed

The same format can also specify CNV regions to the FreeBayes variant caller with FreeBayes’s --cnv-map option:

# Show only CNV regions
cnvkit.py export bed Sample.call.cns -o all-samples.cnv-map.bed

By default only regions with copy number different from the given ploidy (default 2) are output. (Notice what this means for allosomes.) To output all segments, use the --show all option.

vcf¶

Convert segments, ideally already adjusted by the rescale and/or call commands, to a VCF file. Copy ratios are converted to absolute integers, as with BED export, and VCF records are created for the segments where the copy number is different from the expected ploidy (e.g. 2 on autosomes, 1 on haploid sex chromosomes, depending on sample gender).

Gender can be specified with the -g/--gender option, or will be guessed automatically. If a male reference is used, use -y/--male-reference to say so. Note that these are different: If a female sample is run with a male reference, segments on chromosome X with log2-ratio +1 will be skipped, because that’s the expected copy number, while an X-chromosome segment with log2-ratio 0 will be printed as a hemizygous loss.

cnvkit.py export vcf Sample.cns -y -g female -i "SampleID" -o Sample.cnv.vcf

cdt, jtv¶

A collection of probe-level copy ratio files (*.cnr) can be exported to Java TreeView via the standard CDT format or a plain text table:

cnvkit.py export jtv *.cnr -o Samples-JTV.txt
cnvkit.py export cdt *.cnr -o Samples.cdt

seg¶

Similarly, the segmentation files for multiple samples (*.cns) can be exported to the standard SEG format to be loaded in the Integrative Genomic Viewer (IGV):

cnvkit.py export seg *.cns -o Samples.seg

nexus-basic¶

The format nexus-basic can be loaded directly by the commercial program Biodiscovery Nexus Copy Number, specifying the “basic” input format in that program. This allows viewing CNVkit data as if it were from array CGH.

This is a tabular format very similar to .cnr files, with the columns:

1. chromosome
2. start
3. end
4. log2

nexus-ogt¶

The format nexus-ogt can be loaded directly by the commercial program Biodiscovery Nexus Copy Number, specifying the “Custom-OGT” input format in that program. This allows viewing CNVkit data as if it were from a SNP array.

This is a tabular format similar to .cnr files, but with B-allele frequencies (BAFs) extracted from a corresponding VCF file. The format’s columns are (with .cnr equivalents):

1. “Chromosome” (chromosome)
2. “Position” (start)
3. “Position” (end)
4. “Log R Ratio” (log2)
5. “B-Allele Frequency” (from VCF)

The positions of each heterozygous variant record in the given VCF are matched to bins in the given .cnr file, and the variant allele frequencies are extracted and assigned to the matching bins.

• If a bin contains no variants, the BAF field is left blank
• If a bin contains multiple variants, the BAFs of those variants are “mirrored” to be all above .5 (e.g. BAF of .3 becomes .7), then the median is taken as the bin-wide BAF.

Using CNVkit with THetA2¶

CNVkit provides wrappers for exporting segments to THetA2’s input format and importing THetA2’s result file as CNVkit’s segmented .cns files.

THetA2’s input file is a BED-like file, typically with the extension .input, listing the read counts within each copy-number segment in a pair of tumor and normal samples. CNVkit can generate this file given the CNVkit-inferred tumor segmentation (.cns) and normal copy log2-ratios (.cnr) or copy number reference file (.cnn). This bypasses the initial step of THetA2, CreateExomeInput, which counts the reads in each sample’s BAM file.

After running the CNVkit Copy number calling pipeline on a sample, create the THetA2 input file:

# From a paired normal sample
cnvkit.py export theta Sample_Tumor.cns Sample_Normal.cnr -o Sample.theta2.input
# From an existing CNVkit reference
cnvkit.py export theta Sample_Tumor.cns reference.cnn -o Sample.theta2.input

Then, run THetA2 (assuming the program was unpacked at /path/to/theta2/):

# Generates Sample.theta2.BEST.results:
/path/to/theta2/bin/RunTHetA Sample.theta2.input
# Parameters for low-quality samples:
/path/to/theta2/python/RunTHetA.py Sample.theta2.input -n 2 -k 4 -m .90 --FORCE --NUM_PROCESSES `nproc`

Finally, import THetA2’s results back into CNVkit’s .cns format, matching the original segmentation (.cns) to the THetA2-inferred absolute copy number values.:

cnvkit.py import-theta Sample_Tumor.cns Sample.theta2.BEST.results

THetA2 adjusts the segment log2 values to the inferred cellularity of each detected subclone; this can result in one or two .cns files representing subclones if more than one clonal tumor cell population was detected. THetA2 also performs some significance testing of each segment representing a CNA, so there may be fewer segments derived from THetA2 than were originally found by CNVkit.

The segment values are still log2-transformed in the resulting .cns files, for convenience in plotting etc. with CNVkit. These files are also easily converted to other formats using the export command.

gary¶

A generic array of genomic positions.

class cnvlib.gary.GenomicArray(data_table, meta_dict=None)[source]¶

Bases: object

An array of genomic intervals. Base class for genomic data structures.

Can represent most BED-like tabular formats with arbitrary additional columns.

add(other)[source]¶

Combine this array’s data with another GenomicArray (in-place).

Any optional columns must match between both arrays.

add_columns(**columns)[source]¶

Create a new CNA, adding the specified extra columns to this CNA.

as_columns(**columns)[source]¶

Wrap the named columns in this instance’s metadata.

as_dataframe(dframe)[source]¶

Wrap the given pandas dataframe in this instance’s metadata.

as_rows(rows)[source]¶

Wrap the given rows in this instance’s metadata.

autosomes(also=())[source]¶

Select chromosomes w/ integer names, ignoring any ‘chr’ prefixes.

by_chromosome()[source]¶

Iterate over bins grouped by chromosome name.

by_ranges(other, mode='inner', keep_empty=True)[source]¶

Group rows by another GenomicArray’s bin coordinate ranges.

Returns an iterable of (bin, GenomicArray of overlapping rows)). Usually used for grouping probes or SNVs by CNV segments.

mode determines what to do with bins that overlap a boundary of the selection. Values are:

• inner: Drop the bins on the selection boundary, don’t emit them.
• outer: Keep/emit those bins as they are.
• trim: Emit those bins but alter their boundaries to match the selection; the bin start or end position is replaced with the selection boundary position. [default]

Bins in this array that fall outside the other array’s bins are skipped.

chromosome¶

concat(others)[source]¶

Concatenate several GenomicArrays, keeping this array’s metadata.

This array’s data table is not implicitly included in the result.

coords(also=())[source]¶

Iterate over plain coordinates of each bin: chromosome, start, end.

With also, also include those columns.

Example, yielding rows in BED format:

>>> probes.coords(also=["name", "strand"])

copy()[source]¶

Create an independent copy of this object.

drop_extra_columns()[source]¶

Remove any optional columns from this GenomicArray.

Returns a new copy with only the core columns retained:

log2 value, chromosome, start, end, bin name.

end¶

classmethod from_columns(columns, meta_dict=None)[source]¶

Create a new instance from column arrays, given as a dict.

classmethod from_rows(rows, columns=None, meta_dict=None)[source]¶

Create a new instance from a list of rows, as tuples or arrays.

in_range(chrom=None, start=None, end=None, mode='inner')[source]¶

Get the GenomicArray portion within the given genomic range.

mode works as in by_ranges: outer includes bins straddling the range boundaries, trim additionally alters the straddling bins’ endpoints to match the range boundaries, and inner excludes those bins.

in_ranges(chrom=None, starts=None, ends=None, mode='inner')[source]¶

Get the GenomicArray portion within the specified ranges.

Same as in_ranges but the starts and ends are arrays of equal length, and the output concatenates all the selected bins.

keep_columns(columns)[source]¶

Extract a subset of columns, reusing this instance’s metadata.

labels()[source]¶

match_to_bins(other, key, default=0.0, fill=False, summary_func=<function median>)[source]¶

Take values of the other array at each of this array’s bins.

Assign default to indices that fall outside the other array’s bins, or chromosomes that appear in self but not other.

Return an array of the key column values in other corresponding to this array’s bin locations, the same length as this array.

classmethod read(infile, sample_id=None)[source]¶

static row2label(row)[source]¶

sample_id¶

select(selector=None, **kwargs)[source]¶

Take a subset of rows where the given condition is true.

Arguments can be a function (lambda expression) returning a bool, which will be used to select True rows, and/or keyword arguments like gene=”Background” or chromosome=”chr7”, which will select rows where the keyed field equals the specified value.

shuffle()[source]¶

Randomize the order of bins in this array (in-place).

sort()[source]¶

Sort this array’s bins in-place, with smart chromosome ordering.

sort_columns()[source]¶

Sort this array’s columns in-place, per class definition.

start¶

write(outfile=None)[source]¶

Write the wrapped data table to a file or handle in tabular format.

The format is BED-like, but with a header row included and with arbitrary extra columns.

To combine multiple samples in one file and/or convert to another format, see the ‘export’ subcommand.

cnary¶

CNVkit’s core data structure, a copy number array.

class cnvlib.cnary.CopyNumArray(data_table, meta_dict=None)[source]¶

Bases: cnvlib.gary.GenomicArray

An array of genomic intervals, treated like aCGH probes.

Required columns: chromosome, start, end, gene, log2

Optional columns: gc, rmask, spread, weight, probes

by_gene(ignore=('-', '.', 'CGH'))[source]¶

Iterate over probes grouped by gene name.

Emits pairs of (gene name, CNA of rows with same name)

Groups each series of intergenic bins as a ‘Background’ gene; any ‘Background’ bins within a gene are grouped with that gene. Bins with names in ignore are treated as ‘Background’ bins, but retain their name.

center_all(estimator=<function median>, skip_low=False)[source]¶

Recenter coverage values to the autosomes’ average (in-place).

drop_low_coverage()[source]¶

Drop bins with extremely low log2 coverage values.

These are generally bins that had no reads mapped, and so were substituted with a small dummy log2 value to avoid divide-by-zero errors.

expect_flat_cvg(is_male_reference=None)[source]¶

Get the uninformed expected copy ratios of each bin.

Create an array of log2 coverages like a “flat” reference.

This is a neutral copy ratio at each autosome (log2 = 0.0) and sex chromosomes based on whether the reference is male (XX or XY).

get_relative_chrx_cvg()[source]¶

Get the relative log-coverage of chrX in a sample.

guess_average_depth(segments=None, window=100)[source]¶

Estimate the effective average read depth from variance.

Assume read depths are Poisson distributed, converting log2 values to absolute counts. Then the mean depth equals the variance , and the average read depth is the estimated mean divided by the estimated variance. Use robust estimators (Tukey’s biweight location and midvariance) to compensate for outliers and overdispersion.

With segments, take the residuals of this array’s log2 values from those of the segments to remove the confounding effect of real CNVs.

If window is an integer, calculate and subtract a smoothed trendline to remove the effect of CNVs without segmentation (skipped if segments are given).

See: http://www.evanmiller.org/how-to-read-an-unlabeled-sales-chart.html

guess_xx(male_reference=False, verbose=True)[source]¶

Guess whether a sample is female from chrX relative coverages.

Recommended cutoff values:

-0.5 – raw target data, not yet corrected +0.5 – probe data already corrected on a male profile

log2¶

residuals(segments=None)[source]¶

Difference in log2 value of each bin from its segment mean.

If segments are just regions (e.g. RegionArray) with no log2 values precalculated, subtract the median of this array’s log2 values within each region. If no segments are given, subtract each chromosome’s median.

shift_xx(male_reference=False)[source]¶

Adjust chrX coverages (divide in half) for apparent female samples.

squash_genes(summary_func=<function biweight_location>, squash_background=False, ignore=('-', '.', 'CGH'))[source]¶

Combine consecutive bins with the same targeted gene name.

The ignore parameter lists bin names that not be counted as genes to be output.

Parameter summary_func is a function that summarizes an array of coverage values to produce the “squashed” gene’s coverage value. By default this is the biweight location, but you might want median, mean, max, min or something else in some cases.

rary¶

An array of genomic regions or features.

class cnvlib.rary.RegionArray(data_table, meta_dict=None)[source]¶

Bases: cnvlib.gary.GenomicArray

An array of genomic intervals.

classmethod read(fname, sample_id=None, fmt=None)[source]¶

Read regions in any of the expected file formats.

Iterates over tuples of the tabular contents. Header lines are skipped.

Start and end coordinates are base-0, half-open.

write(outfile=<open file '<stdout>', mode 'w'>, fmt='bed', verbose=True)[source]¶

vary¶

An array of genomic intervals, treated as variant loci.

class cnvlib.vary.VariantArray(data_table, meta_dict=None)[source]¶

Bases: cnvlib.gary.GenomicArray

An array of genomic intervals, treated as variant loci.

Required columns: chromosome, start, end, ref, alt

mirrored_baf(above_half=True, tumor_boost=True)[source]¶

Mirrored B-allele frequencies (BAFs).

Flip BAFs to be all above 0.5 (if above_half) or below 0.5, for consistency.

With tumor_boost, normalize tumor-sample allele frequencies to the matched normal allele frequencies.

classmethod read_vcf(infile, sample_id=None, normal_id=None, min_depth=None, skip_hom=False, skip_reject=False, skip_somatic=False)[source]¶

Parse SNV coordinates from a VCF file into a VariantArray.

tumor_boost()[source]¶

TumorBoost normalization of tumor-sample allele frequencies.

De-noises the signal for detecting LOH.

commands¶

The public API for each of the commands defined in the CNVkit workflow.

Command-line interface and corresponding API for CNVkit.

cnvlib.commands.batch_make_reference(normal_bams, target_bed, antitarget_bed, male_reference, fasta, annotate, short_names, split, target_avg_size, access, antitarget_avg_size, antitarget_min_size, output_reference, output_dir, processes, by_count)[source]¶

Build the CN reference from normal samples, targets and antitargets.

cnvlib.commands.batch_run_sample(bam_fname, target_bed, antitarget_bed, ref_fname, output_dir, male_reference=False, scatter=False, diagram=False, rlibpath=None, by_count=False)[source]¶

Run the pipeline on one BAM file.

cnvlib.commands.batch_write_coverage(bed_fname, bam_fname, out_fname, by_count)[source]¶

Run coverage on one sample, write to file.

cnvlib.commands.csvstring(text)[source]¶

cnvlib.commands.do_access(fa_fname, exclude_fnames=(), min_gap_size=5000)[source]¶

List the locations of accessible sequence regions in a FASTA file.

cnvlib.commands.do_antitarget(target_bed, access_bed=None, avg_bin_size=100000, min_bin_size=None)[source]¶

Derive a background/antitarget BED file from a target BED file.

cnvlib.commands.do_breaks(probes, segments, min_probes=1)[source]¶

List the targeted genes in which a copy number breakpoint occurs.

cnvlib.commands.do_call(cnarr, variants=None, method='threshold', ploidy=2, purity=None, is_reference_male=False, is_sample_female=False, thresholds=(-1.1, -0.25, 0.2, 0.7))[source]¶

cnvlib.commands.do_coverage(bed_fname, bam_fname, by_count=False, min_mapq=0)[source]¶

Calculate coverage in the given regions from BAM read depths.

cnvlib.commands.do_fix(target_raw, antitarget_raw, reference, do_gc=True, do_edge=True, do_rmask=True)[source]¶

Combine target and antitarget coverages and correct for biases.

cnvlib.commands.do_gainloss(probes, segments=None, male_reference=False, threshold=0.2, min_probes=3, skip_low=False)[source]¶

Identify targeted genes with copy number gain or loss.

cnvlib.commands.do_heatmap(cnarrs, show_range=None, do_desaturate=False)[source]¶

Plot copy number for multiple samples as a heatmap.

cnvlib.commands.do_import_theta(segarr, theta_results_fname, ploidy=2)[source]¶

cnvlib.commands.do_reference(target_fnames, antitarget_fnames, fa_fname=None, male_reference=False, do_gc=True, do_edge=True, do_rmask=True)[source]¶

Compile a coverage reference from the given files (normal samples).

cnvlib.commands.do_reference_flat(targets, antitargets, fa_fname=None, male_reference=False)[source]¶

Compile a neutral-coverage reference from the given intervals.

Combines the intervals, shifts chrX values if requested, and calculates GC and RepeatMasker content from the genome FASTA sequence.

cnvlib.commands.do_rescale(cnarr, ploidy=2, purity=None, is_reference_male=False, is_sample_female=False)[source]¶

cnvlib.commands.do_scatter(cnarr, segments=None, variants=None, show_range=None, show_gene=None, background_marker=None, do_trend=False, window_width=1000000.0, y_min=None, y_max=None, title=None)[source]¶

Plot probe log2 coverages and CBS calls together.

show_gene: name of gene to highligh show_range: chromosome name or coordinate string like “chr1:20-30”

cnvlib.commands.do_targets(bed_fname, annotate=None, do_short_names=False, do_split=False, avg_size=266.6666666666667)[source]¶

Transform bait intervals into targets more suitable for CNVkit.

cnvlib.commands.parse_args(args=None)[source]¶

Parse the command line.

cnvlib.commands.print_version(_args)[source]¶

Display this program’s version.

cnvlib.commands.rescale_baf(purity, observed_baf, normal_baf=0.5)[source]¶

Adjust B-allele frequencies for sample purity.

Math:

t_baf*purity + n_baf*(1-purity) = obs_baf obs_baf - n_baf * (1-purity) = t_baf * purity t_baf = (obs_baf - n_baf * (1-purity))/purity

cnvlib.commands.verify_gender_arg(cnarr, gender_arg, is_male_reference)[source]¶

The following modules implement lower-level functionality specific to each of the CNVkit sub-commands.

antitarget¶

Supporting functions for the ‘antitarget’ command.

cnvlib.antitarget.find_background_regions(access_chroms, target_chroms, pad_size)[source]¶

Take coordinates of accessible regions and targets; emit antitargets.

cnvlib.antitarget.get_background(target_bed, access_bed, avg_bin_size, min_bin_size)[source]¶

Generate background intervals from target intervals.

Procedure:

• Invert target intervals

• Subtract the inverted targets from accessible regions

• For each of the resulting regions:

  • Shrink by a fixed margin on each end
  • If it’s smaller than min_bin_size, skip
  • Divide into equal-size (region_size/avg_bin_size) portions
  • Emit the (chrom, start, end) coords of each portion

cnvlib.antitarget.guess_chromosome_regions(target_chroms, telomere_size)[source]¶

Determine (minimum) chromosome lengths from target coordinates.

call¶

Call copy number variants from segmented log2 ratios.

cnvlib.call.absolute_clonal(cnarr, ploidy, purity, is_reference_male, is_sample_female)[source]¶

Calculate absolute copy number values from segment or bin log2 ratios.

cnvlib.call.absolute_dataframe(cnarr, ploidy, purity, is_reference_male, is_sample_female)[source]¶

Absolute, expected and reference copy number in a DataFrame.

cnvlib.call.absolute_expect(cnarr, ploidy, is_sample_female)[source]¶

Absolute integer number of expected copies in each bin.

I.e. the given ploidy for autosomes, and XY or XX sex chromosome counts according to the sample’s specified gender.

cnvlib.call.absolute_pure(cnarr, ploidy, is_reference_male)[source]¶

Calculate absolute copy number values from segment or bin log2 ratios.

cnvlib.call.absolute_reference(cnarr, ploidy, is_reference_male)[source]¶

Absolute integer number of reference copies in each bin.

I.e. the given ploidy for autosomes, 1 or 2 X according to the reference gender, and always 1 copy of Y.

cnvlib.call.absolute_threshold(cnarr, ploidy, thresholds, is_reference_male)[source]¶

Call integer copy number using hard thresholds for each level.

Integer values are assigned for log2 ratio values less than each given threshold value in sequence, counting up from zero. Above the last threshold value, integer copy numbers are called assuming full purity, diploidy, and rounding up.

Default thresholds follow this heuristic for calling CNAs in a tumor sample: For single-copy gains and losses, assume 50% tumor cell clonality (including normal cell contamination). Then:

R> log2(2:6 / 4)
-1.0  -0.4150375  0.0  0.3219281  0.5849625

Allowing for random noise of +/- 0.1, the cutoffs are:

DEL(0)  <  -1.1
LOSS(1) <  -0.25
GAIN(3) >=  +0.2
AMP(4)  >=  +0.7

For germline samples, better precision could be achieved with:

LOSS(1) <  -0.4
GAIN(3) >=  +0.3

cnvlib.call.log2_ratios(cnarr, absolutes, ploidy, is_reference_male, min_abs_val=0.001, round_to_int=False)[source]¶

Convert absolute copy numbers to log2 ratios.

Optionally round copy numbers to integers.

Account for reference gender & ploidy of sex chromosomes.

coverage¶

Supporting functions for the ‘antitarget’ command.

cnvlib.coverage.bam_total_reads(bam_fname)[source]¶

Count the total number of mapped reads in a BAM file.

Uses the BAM index to do this quickly.

cnvlib.coverage.bedcov(bed_fname, bam_fname, min_mapq)[source]¶

Calculate depth of all regions in a BED file via samtools (pysam) bedcov.

i.e. mean pileup depth across each region.

cnvlib.coverage.interval_coverages(bed_fname, bam_fname, by_count, min_mapq)[source]¶

Calculate log2 coverages in the BAM file at each interval.

cnvlib.coverage.interval_coverages_count(bed_fname, bam_fname, min_mapq)[source]¶

Calculate log2 coverages in the BAM file at each interval.

cnvlib.coverage.interval_coverages_pileup(bed_fname, bam_fname, min_mapq)[source]¶

Calculate log2 coverages in the BAM file at each interval.

cnvlib.coverage.region_depth_count(bamfile, chrom, start, end, min_mapq)[source]¶

Calculate depth of a region via pysam count.

i.e. counting the number of read starts in a region, then scaling for read length and region width to estimate depth.

Coordinates are 0-based, per pysam.

diagram¶

Chromosome diagram drawing functions.

This uses and abuses Biopython’s BasicChromosome module. It depends on ReportLab, too, so we isolate this functionality here so that the rest of CNVkit will run without it. (And also to keep the codebase tidy.)

cnvlib.diagram.bc_chromosome_draw_label(self, cur_drawing, label_name)[source]¶

Monkeypatch to Bio.Graphics.BasicChromosome.Chromosome._draw_label.

Draw a label for the chromosome. Mod: above the chromosome, not below.

cnvlib.diagram.bc_organism_draw(org, title, wrap=12)[source]¶

Modified copy of Bio.Graphics.BasicChromosome.Organism.draw.

Instead of stacking chromosomes horizontally (along the x-axis), stack rows vertically, then proceed with the chromosomes within each row.

Arguments:

• title: The output title of the produced document.

cnvlib.diagram.build_chrom_diagram(features, chr_sizes, sample_id)[source]¶

Create a PDF of color-coded features on chromosomes.

cnvlib.diagram.create_diagram(cnarr, segarr, threshold, min_probes, outfname, male_reference)[source]¶

Create the diagram.

export¶

Export CNVkit objects and files to other formats.

cnvlib.export.create_chrom_ids(segments)[source]¶

Map chromosome names to integers in the order encountered.

cnvlib.export.export_bed(segments, ploidy, is_reference_male, is_sample_female, label, show)[source]¶

Convert a copy number array to a BED-like DataFrame.

For each region in each sample (possibly filtered according to show), the columns are:

• reference sequence name
• start (0-indexed)
• end
• sample name or given label
• integer copy number

By default (show=”ploidy”), skip regions where copy number is the default ploidy, i.e. equal to 2 or the value set by –ploidy. If show=”variant”, skip regions where copy number is neutral, i.e. equal to the reference ploidy on autosomes, or half that on sex chromosomes.

cnvlib.export.export_nexus_basic(sample_fname)[source]¶

Biodiscovery Nexus Copy Number “basic” format.

Only represents one sample per file.

cnvlib.export.export_nexus_ogt(sample_fname, vcf_fname)[source]¶

Biodiscovery Nexus Copy Number “Custom-OGT” format.

To create the b-allele frequencies column, alterate allele frequencies from the VCF are aligned to the .cnr file bins. Bins that contain no variants are left blank; if a bin contains multiple variants, then the frequencies are all “mirrored” to be above .5, then the median of those values is taken.

cnvlib.export.export_seg(sample_fnames)[source]¶

SEG format for copy number segments.

Segment breakpoints are not the same across samples, so samples are listed in serial with the sample ID as the left column.

cnvlib.export.export_theta(tumor, reference)[source]¶

Convert tumor segments and normal .cnr or reference .cnn to THetA input.

Follows the THetA segmentation import script but avoid repeating the pileups, since we already have the mean depth of coverage in each target bin.

The options for average depth of coverage and read length do not matter crucially for proper operation of THetA; increased read counts per bin simply increase the confidence of THetA’s results.

THetA2 input format is tabular, with columns:

ID, chrm, start, end, tumorCount, normalCount

where chromosome IDs (“chrm”) are integers 1 through 24.

cnvlib.export.export_vcf(segments, ploidy, is_reference_male, is_sample_female, sample_id=None)[source]¶

Convert segments to Variant Call Format.

For now, only 1 sample per VCF. (Overlapping CNVs seem tricky.)

Spec: https://samtools.github.io/hts-specs/VCFv4.2.pdf

cnvlib.export.fmt_cdt(sample_ids, table)[source]¶

Format as CDT.

cnvlib.export.fmt_gct(sample_ids, table)[source]¶

cnvlib.export.fmt_jtv(sample_ids, table)[source]¶

Format for Java TreeView.

cnvlib.export.format_theta_row(seg, cnarr, chrom_id, log2_to_count)[source]¶

Convert a segment’s info to a row of THetA input.

For the normal/reference bin count, take the mean of the bin values within each segment so that segments match between tumor and normal.

cnvlib.export.merge_samples(filenames)[source]¶

Merge probe values from multiple samples into a 2D table (of sorts).

Input:

dict of {sample ID: (probes, values)}

Output:

list-of-tuples: (probe, log2 coverages...)

cnvlib.export.mirrored_baf_median(vals)[source]¶

cnvlib.export.segments2vcf(segments, ploidy, is_reference_male, is_sample_female)[source]¶

Convert copy number segments to VCF records.

fix¶

Supporting functions for the ‘fix’ command.

cnvlib.fix.apply_weights(cnarr, ref_matched, epsilon=0.0001)[source]¶

Calculate weights for each bin.

Weights are derived from:

• bin sizes
• average bin coverage depths in the reference
• the “spread” column of the reference.

cnvlib.fix.center_by_window(cnarr, fraction, sort_key)[source]¶

Smooth out biases according to the trait specified by sort_key.

E.g. correct GC-biased probes by windowed averaging across similar-GC probes; or for similar interval sizes.

cnvlib.fix.edge_gains(target_sizes, gap_sizes, insert_size)[source]¶

Calculate coverage gain from neighboring baits’ flanking reads.

Letting i = insert size, t = target size, g = gap to neighboring bait, the gain of coverage due to a nearby bait, if g < i, is:

(i-g)^2 / 4it

If the neighbor flank extends beyond the target (t+g < i), reduce by:

(i-t-g)^2 / 4it

If a neighbor overlaps the target, treat it as adjacent (gap size 0).

cnvlib.fix.edge_losses(target_sizes, insert_size)[source]¶

Calculate coverage losses at the edges of baited regions.

Letting i = insert size and t = target size, the proportional loss of coverage near the two edges of the baited region (combined) is:

i/2t

If the “shoulders” extend outside the bait $(t < i), reduce by:

(i-t)^2 / 4it

on each side, or (i-t)^2 / 2it total.

cnvlib.fix.get_edge_bias(cnarr, margin)[source]¶

Quantify the “edge effect” of the target tile and its neighbors.

The result is proportional to the change in the target’s coverage due to these edge effects, i.e. the expected loss of coverage near the target edges and, if there are close neighboring tiles, gain of coverage due to “spill over” reads from the neighbor tiles.

(This is not the actual change in coverage. This is just a tribute.)

cnvlib.fix.load_adjust_coverages(pset, ref_pset, skip_low, fix_gc, fix_edge, fix_rmask)[source]¶

Load and filter probe coverages; correct using reference and GC.

cnvlib.fix.mask_bad_probes(probes)[source]¶

Flag the probes with excessively low or inconsistent coverage.

Returns a bool array where True indicates probes that failed the checks.

cnvlib.fix.match_ref_to_probes(ref_pset, probes)[source]¶

Filter the reference probes to match the target or antitarget probe set.

importers¶

Import from other formats to the CNVkit format.

cnvlib.importers.find_picard_files(file_and_dir_names)[source]¶

Search the given paths for ‘targetcoverage’ CSV files.

Per the convention we use in our Picard applets, the target coverage file names end with ‘.targetcoverage.csv’; anti-target coverages end with ‘.antitargetcoverage.csv’.

cnvlib.importers.import_picard_pertargetcoverage(fname)[source]¶

Parse a Picard CalculateHsMetrics PER_TARGET_COVERAGE file.

Return a CopyNumArray.

Input column names:

chrom (str), start, end, length (int), name (str), %gc, mean_coverage, normalized_coverage (float)

cnvlib.importers.import_seg(segfname, chrom_names, chrom_prefix, from_log10)[source]¶

Parse a SEG file as an iterable of CopyNumArray instances.

chrom_names:

Map (string) chromosome IDs to names. (Applied before chrom_prefix.) e.g. {‘23’: ‘X’, ‘24’: ‘Y’, ‘25’: ‘M’}

chrom_prefix: prepend this string to chromosome names

(usually ‘chr’ or None)

from_log10: Convert values from log10 to log2.

cnvlib.importers.parse_theta_results(fname)[source]¶

Parse THetA results into a data structure.

Columns: NLL, mu, C, p*

cnvlib.importers.unpipe_name(name)[source]¶

Fix the duplicated gene names Picard spits out.

Return a string containing the single gene name, sans duplications and pipe characters.

Picard CalculateHsMetrics combines the labels of overlapping intervals by joining all labels with ‘|’, e.g. ‘BRAF|BRAF’ – no two distinct targeted genes actually overlap, though, so these dupes are redundant. Meaningless target names are dropped, e.g. ‘CGH|FOO|-‘ resolves as ‘FOO’. In case of ambiguity, the longest name is taken, e.g. “TERT|TERT Promoter” resolves as “TERT Promoter”.

reference¶

Supporting functions for the ‘reference’ command.

cnvlib.reference.bed2probes(bed_fname)[source]¶

Create neutral-coverage probes from intervals.

cnvlib.reference.calculate_gc_lo(subseq)[source]¶

Calculate the GC and lowercase (RepeatMasked) content of a string.

cnvlib.reference.combine_probes(filenames, fa_fname, is_male_reference, skip_low, fix_gc, fix_edge, fix_rmask)[source]¶

Calculate the median coverage of each bin across multiple samples.

Input:

List of .cnn files, as generated by ‘coverage’ or ‘import-picard’. fa_fname: fil columns for GC and RepeatMasker genomic values.

Returns:

A single CopyNumArray summarizing the coverages of the input samples, including each bin’s “average” coverage, “spread” of coverages, and genomic GC content.

cnvlib.reference.get_fasta_stats(probes, fa_fname)[source]¶

Calculate GC and RepeatMasker content of each bin in the FASTA genome.

cnvlib.reference.reference2regions(reference, coord_only=False)[source]¶

Extract iterables of target and antitarget regions from a reference.

Like loading two BED files with ngfrills.parse_regions.

cnvlib.reference.warn_bad_probes(probes)[source]¶

Warn about target probes where coverage is poor.

Prints a formatted table to stderr.

reports¶

Supports the sub-commands breaks and gainloss.

Supporting functions for the text/tabular-reporting commands.

Namely: breaks, gainloss.

cnvlib.reports.gainloss_by_gene(probes, threshold, skip_low=False)[source]¶

Identify genes where average bin copy ratio value exceeds threshold.

NB: Must shift sex-chromosome values beforehand with shift_xx, otherwise all chrX/chrY genes may be reported gained/lost.

cnvlib.reports.gainloss_by_segment(probes, segments, threshold, skip_low=False)[source]¶

Identify genes where segmented copy ratio exceeds threshold.

NB: Must shift sex-chromosome values beforehand with shift_xx, otherwise all chrX/chrY genes may be reported gained/lost.

cnvlib.reports.get_breakpoints(intervals, segments, min_probes)[source]¶

Identify CBS segment breaks within the targeted intervals.

cnvlib.reports.get_gene_intervals(all_probes, ignore=('-', '.', 'CGH'))[source]¶

Tally genomic locations of each targeted gene.

Return a dict of chromosomes to a list of tuples: (gene name, start, end).

cnvlib.reports.group_by_genes(probes, skip_low)[source]¶

Group probe and coverage data by gene.

Return an iterable of genes, in chromosomal order, associated with their location and coverages:

[(gene, chrom, start, end, [coverages]), ...]

segmentation¶

Segmentation of copy number values.

cnvlib.segmentation.do_segmentation(cnarr, method, threshold=None, variants=None, skip_low=False, skip_outliers=10, save_dataframe=False, rlibpath=None)[source]¶

Infer copy number segments from the given coverage table.

cnvlib.segmentation.drop_outliers(cnarr, width, factor)[source]¶

Drop outlier bins with log2 ratios too far from the trend line.

Outliers are the log2 values factor times the 90th quantile of absolute deviations from the rolling average, within a window of given width. The 90th quantile is about 1.97 standard deviations if the log2 values are Gaussian, so this is similar to calling outliers factor * 1.97 standard deviations from the rolling mean. For a window size of 50, the breakdown point is 2.5 outliers within a window, which is plenty robust for our needs.

cnvlib.segmentation.repair_segments(segments, orig_probes)[source]¶

Post-process segmentation output.

1. Ensure every chromosome has at least one segment.
2. Ensure first and last segment ends match 1st/last bin ends (but keep log2 as-is).
3. Store probe-level gene names, comma-separated, as the segment name.

cnvlib.segmentation.seg2cns(seg_text)[source]¶

Convert R dataframe contents (SEG) to our native tabular format.

Return a pandas.Dataframe with CNA columns.

cnvlib.segmentation.squash_segments(seg_pset)[source]¶

Combine contiguous segments.

cnvlib.segmentation.transfer_names_weights(segments, cnarr, ignore=('-', '.', 'CGH'))[source]¶

Copy gene names from cnarr to the segmented segarr.

Segment name is the comma-separated list of bin gene names.

target¶

Transform bait intervals into targets more suitable for CNVkit.

cnvlib.target.add_refflat_names(region_rows, refflat_fname)[source]¶

Apply RefSeq gene names to a list of targeted regions.

cnvlib.target.assign_names(region_rows, refflat_fname, default_name='-')[source]¶

Replace the interval gene names with those at the same loc in refFlat.txt

cnvlib.target.emit(coords, names)[source]¶

Try filtering names again. Format the row for yielding.

cnvlib.target.filter_names(names, exclude=('mRNA', ))[source]¶

Remove less-meaningful accessions from the given set.

cnvlib.target.parse_refflat_line(line)[source]¶

Parse one line of refFlat.txt; return relevant info.

Pair up the exon start and end positions. Add strand direction to each chromosome as a “+”/”-” suffix (it will be stripped off later).

cnvlib.target.read_refflat_genes(fname)[source]¶

Parse genes; merge those with same name and overlapping regions.

Returns a dict of: {(chrom, strand): [(gene start, gene end, gene name), ...]}

cnvlib.target.shorten_labels(interval_rows)[source]¶

Reduce multi-accession interval labels to the minimum consistent.

So: BED or interval_list files have a label for every region. We want this to be a short, unique string, like the gene name. But if an interval list is instead a series of accessions, including additional accessions for sub-regions of the gene, we can extract a single accession that covers the maximum number of consecutive regions that share this accession.

e.g.

... chr1 879125 879593 + mRNA|JX093079,ens|ENST00000342066,mRNA|JX093077,ref|SAMD11,mRNA|AF161376,mRNA|JX093104 chr1 880158 880279 + ens|ENST00000483767,mRNA|AF161376,ccds|CCDS3.1,ref|NOC2L ...

becomes:

chr1 879125 879593 + mRNA|AF161376 chr1 880158 880279 + mRNA|AF161376

cnvlib.target.shortest_name(names)[source]¶

Return the shortest trimmed name from the given set.

cnvlib.target.split_targets(region_rows, avg_size)[source]¶

Split large tiled intervals into smaller, consecutive targets.

For each of the tiled regions:

• Divide into equal-size (tile_size/avg_size) portions
• Emit the (chrom, start, end) coords of each portion

Bin the regions according to avg_size.

core¶

CNV utilities.

cnvlib.core.assert_equal(msg, **values)[source]¶

Evaluate and compare two or more values for equality.

Sugar for a common assertion pattern. Saves re-evaluating (and retyping) the same values for comparison and error reporting.

Example:

>>> assert_equal("Mismatch", expected=1, saw=len(['xx', 'yy']))
...
ValueError: Mismatch: expected = 1, saw = 2

cnvlib.core.check_unique(items, title)[source]¶

Ensure all items in an iterable are identical; return that one item.

cnvlib.core.fbase(fname)[source]¶

Strip directory and all extensions from a filename.

cnvlib.core.rbase(fname)[source]¶

Strip directory and final extension from a filename.

cnvlib.core.sorter_chrom(label)[source]¶

Create a sorting key from chromosome label.

Sort by integers first, then letters or strings. The prefix “chr” (case-insensitive), if present, is stripped automatically for sorting.

E.g. chr1 < chr2 < chr10 < chrX < chrY < chrM

cnvlib.core.sorter_chrom_at(index)[source]¶

Create a sort key function that gets chromosome label at a list index.

cnvlib.core.write_dataframe(outfname, dframe, header=True)[source]¶

Write a pandas.DataFrame to a tabular file.

cnvlib.core.write_text(outfname, text, *more_texts)[source]¶

Write one or more strings (blocks of text) to a file.

cnvlib.core.write_tsv(outfname, rows, colnames=None)[source]¶

Write rows, with optional column header, to tabular file.

metrics¶

Robust estimators of central tendency and scale.

For use in evaluating performance of copy number estimation.

See:

https://en.wikipedia.org/wiki/Robust_measures_of_scale https://astropy.readthedocs.org/en/latest/_modules/astropy/stats/funcs.html

cnvlib.metrics.biweight_location(a, initial=None, c=6.0, epsilon=0.0001)[source]¶

Compute the biweight location for an array.

The biweight is a robust statistic for determining the central location of a distribution.

cnvlib.metrics.biweight_midvariance(a, initial=None, c=9.0, epsilon=0.0001)[source]¶

Compute the biweight midvariance for an array.

The biweight midvariance is a robust statistic for determining the midvariance (i.e. the standard deviation) of a distribution.

See: https://en.wikipedia.org/wiki/Robust_measures_of_scale#The_biweight_midvariance https://astropy.readthedocs.org/en/latest/_modules/astropy/stats/funcs.html

cnvlib.metrics.confidence_interval_bootstrap(bins, alpha, bootstraps=100)[source]¶

Confidence interval for segment mean log2 value, estimated by bootstrap.

cnvlib.metrics.ests_of_scale(deviations)[source]¶

Estimators of scale: standard deviation, MAD, biweight midvariance.

Calculates all of these values for an array of deviations and returns them as a tuple.

cnvlib.metrics.interquartile_range(a)[source]¶

Compute the difference between the array’s first and third quartiles.

cnvlib.metrics.median_absolute_deviation(a, scale_to_sd=True)[source]¶

Compute the median absolute deviation (MAD) of array elements.

The MAD is defined as: median(abs(a - median(a))).

See: https://en.wikipedia.org/wiki/Median_absolute_deviation

cnvlib.metrics.modal_location(arr)[source]¶

Return the modal value of an array’s values.

The “mode” is the location of peak density among the values, estimated using a Gaussian kernel density estimator.

arr is a 1-D array of floating-point values, e.g. bin log2 ratio values.

cnvlib.metrics.prediction_interval(bins, alpha)[source]¶

Prediction interval, estimated by percentiles.

cnvlib.metrics.q_n(a)[source]¶

Rousseeuw & Croux’s (1993) Q_n, an alternative to MAD.

Qn := Cn first quartile of (|x_i - x_j|: i < j)

where Cn is a constant depending on n.

Finite-sample correction factors must be used to calibrate the scale of Qn for small-to-medium-sized samples.

n E[Qn] – —– 10 1.392 20 1.193 40 1.093 60 1.064 80 1.048 100 1.038 200 1.019

cnvlib.metrics.segment_mean(cnarr, skip_low=False)[source]¶

Weighted average of bin log2 values.

ngfrills¶

NGS utilities.

cnvlib.ngfrills.call_quiet(*args)[source]¶

Safely run a command and get stdout; print stderr if there’s an error.

Like subprocess.check_output, but silent in the normal case where the command logs unimportant stuff to stderr. If there is an error, then the full error message(s) is shown in the exception message.

cnvlib.ngfrills.ensure_path(fname)[source]¶

Create dirs and move an existing file to avoid overwriting, if necessary.

If a file already exists at the given path, it is renamed with an integer suffix to clear the way.

cnvlib.ngfrills.safe_write(*args, **kwds)[source]¶

Write to a filename or file-like object with error handling.

If given a file name, open it. If the path includes directories that don’t exist yet, create them. If given a file-like object, just pass it through.

cnvlib.ngfrills.temp_write_text(*args, **kwds)[source]¶

Save text to a temporary file.

NB: This won’t work on Windows b/c the file stays open.

parallel¶

Utilities for multi-core parallel processing.

class cnvlib.parallel.SerialPool[source]¶

Bases: object

Mimic the multiprocessing.Pool interface, but run in serial.

apply_async(func, args)[source]¶

Just call the function.

close()[source]¶

join()[source]¶

cnvlib.parallel.pick_pool(nprocs)[source]¶

params¶

Defines several constants used in the pipeline.

Hard-coded parameters for CNVkit. These should not change between runs.

plots¶

Plotting utilities.

cnvlib.plots.chromosome_sizes(probes, to_mb=False)[source]¶

Create an ordered mapping of chromosome names to sizes.

cnvlib.plots.cnv_on_chromosome(axis, probes, segments, genes, background_marker=None, do_trend=False, y_min=None, y_max=None)[source]¶

Draw a scatter plot of probe values with CBS calls overlaid.

Argument ‘genes’ is a list of tuples: (start, end, gene name)

cnvlib.plots.cnv_on_genome(axis, probes, segments, pad, do_trend=False, y_min=None, y_max=None)[source]¶

Plot coverages and CBS calls for all chromosomes on one plot.

cnvlib.plots.cvg2rgb(cvg, desaturate)[source]¶

Choose a shade of red or blue representing log2-coverage value.

cnvlib.plots.gene_coords_by_name(probes, names)[source]¶

Find the chromosomal position of each named gene in probes.

Returns a dict: {chromosome: [(start, end, gene name), ...]}

cnvlib.plots.gene_coords_by_range(probes, chrom, start, end, ignore=('-', '.', 'CGH'))[source]¶

Find the chromosomal position of all genes in a range.

Returns a dict: {chromosome: [(start, end, gene), ...]}

cnvlib.plots.group_snvs_by_segments(snv_posns, snv_freqs, segments, chrom=None)[source]¶

Group SNP allele frequencies by segment.

Return an iterable of: start, end, value.

cnvlib.plots.parse_range_text(text)[source]¶

Parse a chromosomal range specification.

Range spec string should look like chr1:1234-5678 or chr1:1234- or chr1:-5678, where missing start becomes 0 and missing end becomes None.

cnvlib.plots.partition_by_chrom(chrom_snvs)[source]¶

Group the tumor shift values by chromosome (for statistical testing).

cnvlib.plots.plot_x_dividers(axis, chrom_sizes, pad)[source]¶

Plot vertical dividers and x-axis labels given the chromosome sizes.

Returns a table of the x-position offsets of each chromosome.

Draws vertical black lines between each chromosome, with padding. Labels each chromosome range with the chromosome name, centered in the region, under a tick. Sets the x-axis limits to the covered range.

cnvlib.plots.setup_chromosome(axis, probes=None, segments=None, variants=None, y_min=None, y_max=None, y_label=None)[source]¶

Configure axes for plotting a single chromosome’s data.

probes, segments, and variants should already be subsetted to the region that will be plotted.

cnvlib.plots.setup_genome(axis, probes, segments, variants, y_min=None, y_max=None)[source]¶

Configure axes for plotting a whole genomes’s data.

cnvlib.plots.snv_on_chromosome(axis, variants, segments, genes, do_trend, do_boost=False)[source]¶

cnvlib.plots.snv_on_genome(axis, variants, chrom_sizes, segments, do_trend, pad, do_boost=False)[source]¶

Plot a scatter-plot of SNP chromosomal positions and shifts.

cnvlib.plots.test_loh(bins, alpha=0.0025)[source]¶

Test each chromosome’s SNP shifts and the combined others’.

The statistical test is Mann-Whitney, a one-sided non-parametric test for difference in means.

cnvlib.plots.unpack_range(a_range)[source]¶

Extract chromosome, start, end from a string or tuple.

Examples:

“chr1” -> (“chr1”, None, None) “chr1:100-123” -> (“chr1”, 100, 123) (“chr1”, 100, 123) -> (“chr1”, 100, 123)

smoothing¶

Signal smoothing functions.

cnvlib.smoothing.check_inputs(x, width)[source]¶

Transform width into a half-window size.

width is either a fraction of the length of x or an integer size of the whole window. The output half-window size is truncated to the length of x if needed.

cnvlib.smoothing.fit_edges(x, y, wing, polyorder=3)[source]¶

Apply polynomial interpolation to the edges of y, in-place.

Calculates a polynomial fit (of order polyorder) of x within a window of width twice wing, then updates the smoothed values y in the half of the window closest to the edge.

cnvlib.smoothing.outlier_iqr(a, c=3.0)[source]¶

Detect outliers as a multiple of the IQR from the median.

By convention, “outliers” are points more than 1.5 * IQR from the median, and “extremes” or extreme outliers are those more than 3.0 * IQR.

cnvlib.smoothing.outlier_mad_median(a)[source]¶

MAD-Median rule for detecting outliers.

Returns: a boolean array of the same size, where outlier indices are True.

X_i is an outlier if:

 | X_i - M |
_____________  > K ~= 2.24

 MAD / 0.6745

where $K = sqrt( X^2_{0.975,1} )$, the square root of the 0.975 quantile of a chi-squared distribution with 1 degree of freedom.

This is a very robust rule with the highest possible breakdown point of 0.5.

See:

• Davies & Gather (1993) The Identification of Multiple Outliers.
• Rand R. Wilcox (2012) Introduction to robust estimation and hypothesis testing. Ch.3: Estimating measures of location and scale.

cnvlib.smoothing.rolling_median(x, width)[source]¶

Rolling median with mirrored edges.

cnvlib.smoothing.rolling_outlier_iqr(x, width, c=3.0)[source]¶

Detect outliers as a multiple of the IQR from the median.

By convention, “outliers” are points more than 1.5 * IQR from the median (~2 SD if values are normally distributed), and “extremes” or extreme outliers are those more than 3.0 * IQR (~4 SD).

cnvlib.smoothing.rolling_outlier_quantile(x, width, q, m)[source]¶

Return a boolean mask of outliers by multiples of a quantile in a window.

Outliers are the array elements outside m times the q‘th quantile of deviations from the smoothed trend line, as calculated from the trend line residuals. (For example, take the magnitude of the 95th quantile times 5, and mark any elements greater than that value as outliers.)

This is the smoothing method used in BIC-seq (doi:10.1073/pnas.1110574108) with the parameters width=200, q=.95, m=5 for WGS.

cnvlib.smoothing.rolling_outlier_std(x, width, stdevs)[source]¶

Return a boolean mask of outliers by stdev within a rolling window.

Outliers are the array elements outside stdevs standard deviations from the smoothed trend line, as calculated from the trend line residuals.

cnvlib.smoothing.rolling_quantile(x, width, quantile)[source]¶

Rolling quantile (0–1) with mirrored edges.

cnvlib.smoothing.rolling_std(x, width)[source]¶

Rolling quantile (0–1) with mirrored edges.

cnvlib.smoothing.smoothed(x, width, do_fit_edges=False)[source]¶

Smooth the values in x with the Kaiser windowed filter.

See: https://en.wikipedia.org/wiki/Kaiser_window

Parameters:

x : array-like

1-dimensional numeric data set.

width : float

Fraction of x’s total length to include in the rolling window (i.e. the proportional window width), or the integer size of the window.