CoNVaDING
CoNVaDING User Guide
Written for release CoNVaDING v1.2.0
Introduction
CoNVaDING (Copy Number Variation Detection In Next-generation sequencing Gene panels) was designed for small (single-exon) copy number variation (CNV) detection in high coverage next-generation sequencing (NGS) data, such as obtained by analysis of smaller targeted gene panels.
CoNVaDING makes use of a group of (at least 30) possible control samples from which the samples with the most similar overall pattern are selected as control samples. These samples are then used for read-depth normalization on all (autosomal) targets and on all targets per gene. CNV prediction is based on a combination of ratio scores and Z-scores of the sample of interest compared to the selected controlsamples.
Quality (QC) metrics are calculated per sample and per analyzed target. Output is generated on three levels:
longlist: This list contains all calls, disregarding the target quality.
shortlist: This list contains a subset of the longlist, filtered on within sample target QC metrics.
final list: This list contains a subset of the shortlist, filtered on target QC metrics obtained from other samples.
CoNVaDING has been written for use of CNV detection in high coverage NGS data (at least ~200x). With lower coverages it might still work, but more targets will fail QC metrics.
The program is written in perl and has dependencies on specific perl libraries as well as on samtools version 1.3 or higher.
Installation
The latest version of CoNVaDING can be downloaded from the github page.
CoNVaDING has several dependencies:
The Statistics::Normality package for perl
Samtools (v1.3 or higher)
After installation the samtools executable has to be added to your local environment using the $PATH variable.
This version is known to be compatible with SAMtools version 1.3
General comments
Targets, or the region of interest, usually consist of an exon with some flanking bases, but can be anything that you specify, for instance part of an exon.
Because of the fact that the average coverage of the specified targets is used for the calculations, the resolution of the calls is also dependent on those targets. This means that the exact break point can not be determined. It might be within the last target of the call, or in the intron flanking this exon. Sub-exonic CNV's can also be missed if they if they do not have enough effect on the average coverage of the target.
Analysis
The analysis consists of three steps that have to be run separately and is based on a list of targets in a bed format.
CoNVaDING can be started using the following command:
If no options are used the help menu will appear.
Create normalized count files
The first step in the analysis is to create normalized count files. This can be done in two ways, from a bam file or from a text file including mean coverage per target.
StartWithBam
If a bam file is used CoNVaDING will use samtools to calculate the mean coverage for each target. For this type of analysis the StartWithBam mode has to be selected.
The basic analysis starts as follows:
All bamfiles should be stored in the same folder, which can be specified with the 'inputDir' option The 'outputDir' option should specify the path to the folder in which the normalized coverage files .aligned.only.normalized.coverage.txt should be stored.
The bed file should contain the regions of interest seperated in four columns specifying the chromosome, start position, stop position and the gene. No headers should be included. Make sure the style of the bed file matches the style of the bam file. UCSC style files show 'chr1, chr2, etc.' in the header, while in other styles this is '1, 2, etc.'. The naming of the chromosomes in the bed file should be the same as in the bam file header.
It is important that the gene column has the exact same gene name for every target of the same gene, because these names are used to cluster targets for a normalization based on the targets belonging to the same gene. The bedfile should be sorted on chromosome and start position.
Bed file example:
The analysis options can be further extended:
If a control set is not yet present, or if the samples that are analyzed have to be added to the control set the following options should be added:
The 'useSampleAsControl' option specifies that the samples have to be used as a control sample later on. The 'controlsDir' is the location where the normalized coverage files of the control samples will be stored.
If duplicates have to be removed before coverage calculations use the following option:
This is advisable for capturing data, but should not be done for amplicon data.
StartWithAvgCount
If no bam files are present the analysis can also start with a text file specifying average counts per target. This file shoud contain the headers as shown in the example below. This enables the use of alternative analysis software.
An example of a text file that can be used in this mode:
When this mode is used the analysis is started as follows:
Also here the following options can be used if the samples should be used as a control set in later steps:
The -rmdup option is not available in this mode. If necessary the duplicates should have been removed before calculating the mean coverage per target.
Selecting the most informative control samples
StartWithMatchScore
The next step in the analysis is selecting the control samples estimated to be the most informative.
If a bam file is used CoNVaDING will use samtools to calculate the mean coverage for each target. For this type of analysis the StartWithBam mode has to be selected.
The basic analysis starts as follows:
The 'inputDir' option should specify the path to the folder in which the normalized coverage files .aligned.only.normalized.coverage.txt are be stored (the outputfolder of the previous step).
The 'outputDir' option should specify the path to output folder. The script will produce two types of output files:
.best.match.score.txt shows the matchscore and the paths to the selected control samples.
.normalized.autosomal.coverage.all.controls.txt show the normalized coverage for all possible control samples
The 'controlsDir' option should show the directory in which the control samples that have to be used are stored.
The analysis options can be further extended:
On default only targets located on autosomal chromosomes will be analyzed. If some targets are located on the sex chromosomes the following option should be added:
Note that for this option only samples of the same sex as the sample of interest can be used as possible control samples.
On default 30 samples are selected to create the control group. If you wish to use a different number of control samples this can be indicated with the option:
to select for instance the 40 best matching samples.
CNV Detection
StartWithBestScore
The last step in the analysis is the CNV detection itself.
The basic analysis starts as follows:
The 'inputDir' option should specify the path to the folder in which the .best.match.score.txt files are stored (the outputfolder of the previous step).
The 'outputDir' option should specify the path to output folder. The script will produce four types of output files:
.best.score.log show the used control samples, the sample ratio score and the omitted regions for the sample ratio score calculation.
.best.score.longlist.txt contains all calls, regardless of the target quality
.best.score.shortlist.txt contains the high quality calls based on within sample target QC
.best.score.totallist.txt contains information about all targets (ratio scores, Z-scores, QC)
The analysis options can be further extended:
The sample ratio calculation is based on calculation the coefficient of variation of the normalized targets of the sample of interest. In this calculation highly variable targets are excluded. On default a target is considered highly variable if after transforming the normalized target ratio's of all samples in the possible control group 20 percent or more of the samples have a Z-score outside the -3 to 3 range. This percentage can be altered using the 'regionTreshold' option. For a threshold of 30 percent of the samples for instance the following option can be used:
To alter the ratio thresholds when making a call for a deletion of duplication for a region during the analysis, the ratioCutOffLow and ratioCutOffHigh parameters can be used. To apply a threshold of ratio score below 0.65 for a deletion and above 1.4 for duplication use:
The same thresholds for calling a deletion or duplication can also be applied using the Z-score value cutoff. To call a deletion when the Z-score is below -3 or duplication when the Z-score is above 3 use:
To finetune the variant list one can generate a list of targets which in general are of lower quality in all possible controlsamples and apply this as a filter to generate a final list of high quality calls. This can be done by executing two steps:
GenerateTargetQcList
To generate the list of targets and corresponding quality thresholds run:
For this analysis, the same region threshold, ratio cutoffs and Z-score cutoffs as explained above can be altered using their corresponding parameters.
CreateFinalList
To apply the generated list of sample target QCs to the .best.score.shortlist.txt files execute:
To change the percentage of samples in which a target can be labelled as less reliabe, for example in 20 percent of the samples, use the option:
This produces the following output file:
.finallist.txt contains all final calls, basically a filtered shortlist file.
Test dataset
A test dataset can be downloaded from here. The dataset contains one sample and fifty control samples. The scripts described above are also included in the test dataset folder.
StartWithAvgCount
Analysis starts from the StartWithAvgCount step:
Running this step should create a resultsfolder within the Test_dataset folder containing a normalized coverage file.
StartWithMatchScore
Subsequently the best matching samples can be determined using the StartWithMatchScore option:
The best thirty samples should be selected, as shown below:
StartWithBestScore
Now the CNV detection can be performed, using the StartWithBestScore option:
A StartWithBestScore folder is created in the results folder, containing four files.
Sample1.average.counts.best.score.log
Sample1.average.counts.best.score.longlist.txt
Sample1.average.counts.best.score.shortlist.txt
Sample1.average.counts.best.score.totallist.txt
The log file will show the Sample QC: SAMPLE_CV: 0.0749944561392519
and the Match QC: MEAN_AVERAGE_BEST_MATCHSCORE: 0.0748225246685803
The longlist should contain six calls:
Sample1.average.counts.best.score.longlist.txt
The totallist contains the information of all targets and shows the ratio's and Z-scores and the coefficient of variation of each target of the control set ratio's. If the coefficient of variation of the target (AUTO_VC) is too high (above 0.10) the target QC fails and the target is labelled low quality.
Sample1.average.counts.best.score.totallist.txt
Since the ABCC9 exon in which a duplication was detected has a low quality, this exon is filtered from the shortlist
Sample1.average.counts.best.score.shortlist.txt
GenerateTargetQcList
For final filtering the TargetQCList is made using all control samples:
CreateFinalList
Finally the shortlist is filtered usint the targetQClist:
The call of the titin exon had sufficient quality within the analysed sample. However, the exon performed poorly in a large portion of the control samples. Therefore, the call is filtered from the final list, leaving four calls. Notice that the whole gene deletion of DSC2 and the six exon deletion of DSG2 consist of consecutive targets. It is possible that there is one big deletion, containing both genes. However, CoNVaDING will always treat CNVs in different genes as seperate calls.
Sample1.average.counts.best.score.shortlist.finallist.txt
QC Thresholds
If the following thresholds are exceeded using default settings the CNV calling is less reliable. The target ratio is also used to filter calls for the shortlist. Both QC metrics are used for filtering the final list.
Sample CV: 0.09
Target CV: 0.10
FAQ
Which are the criteria of adding a call to the longlist, shortlist or finallist?
All calls passing threshold criteria in the StartWithBestScore function are added to the longlist, regardless of there quality. These will be all targets showing DEL, DUP, HOM_DEL or HOM_DUP in the totallist. This list provides a high sensitivity, but has a relatively low specificity.
The shortlist takes into account the variability of the target in the selected control samples. If the target CV in the used control samples is higher than 0.10 the target is considered LOW_QUALITY, as shown in the totallist. Calls consisting only of low quality targets are filtered from the longlist and are not present on the shortlist.
The finallist is a further filtering of the shortlist. If one or more targets of a call are not considered low quality in the sample, the call will be on the shortlist. However, if all targets of the call are considered low quality in a higher percentage of samples having a SAMPLE_CV lower than 0.09, than specified in the -regionThreshold option (default 20%), the call will not be present on the finallist. This list has the highest specificity, but CNV's in low quality regions can be missed.
If a target has a coverage of 0 CoNVaDING crashes. For this target the normalized coverage file shows -NaN as region coverage.
This is the result of a bug in CoNVaDING v1.1.6. This should be fixed in version 1.2.0.
A call I expect is present on the shortlist, but not on the finallist.
If the sample containing the CNV is present in the ControlsDir a call can be filtered because of itself. The reason for this is that the sample will be sometimes selected as a control sample. The deviating ratio of reads in the (here control) sample can cause a target to be considered as low quality. This might increase the number of control samples considering the target as low quality to above the threshold, causing the call to be filtered.
This can be prevented by performing a preliminary analysis using all the samples in the ControlsDir as input as well as controls for the StartWithMatchScore and StartWithBestScore option. Samples showing calls in the shortlist should be removed from the ControlsDir. After this the real analysis can be started.
What is the difference between LOW_QUALITY and LOW_QUALITY, FAILED?
Both are low quality targets. However, for the LOW_QUALITY, FAILED targets all targets belonging to the same gene are of low quality. In other words, analysis has failed for the entire gene. In LOW_QUALITY targets at least one of the other targets belonging to the same gene has a good quality.
Can CoNVaDING be used on exome data?
The performance of CoNVaDING is mainly dependent on the number of reads covering a target. If the exome coverage is high enough, CoNVaDING should work well. The optimum coverage is 200x or higher. Lower coverages will result in a higher number of low quality targets. Increasing the target size to include more reads in the target, for instance by adding 100bp flanking regions, may increase the quality.
CoNVaDING can take some time to run when used on large bam files.
The analysis speed of the StartWithBam function can be increased by running several instances of CoNVaDING in parallel. All bam files in the inputDir will be processed in sequence. Creating two inputDirs each containing half of the bam files can half the analysis time in case enough recources are available.
All the targets on the X chromosome are of low quality.
Make sure all control samples are of the same gender.
Changelog v1.2.0
The following bugs have been fixed in the update from version 1.1.6 to v1.2.0:
Major
In version 1.1.6 after running the StartWithBam function analyzing a region having a coverage of 0, the normalized coverage would have a value of -NaN. This issue is fixed in version 1.2.0, by adding the -a flag to the samtools depth function. However, this causes a dependency on samtools version 1.3 or higher.
In version 1.1.6 HOM_DEL and HOM_DUP calls were interpreted as separate events as DEL and DUP calls. This could cause calls to be cut into several pieces and create gaps in the calls because of downstream filtering. Now, HOM_DEL calls are interpreted as part of a DEL call as well as a separate call. For instance if exons 1 and 3 have a DEL call in the totallist and exon 2 a HOM_DEL call. A deletion of exons 1 to 3 will be called as well as a homozygous deletion of exon 2. Similarly HOM_DUP events are considered as part of DUP events as well as a separate call. Note that in case all targets of a call are considered HOM_DEL or HOM_DUP the call will appear twice on the longlist, shortlist and/or finallist.
Minor
In the StartWithMatchScore function the 'Avg abs diff incl sex' was showing the same value as the 'Avg abs diff autosomal'. Now the values are correct.
A bedfile containing UCSC style notation 'chr1, chr2, etc.' is now be compatible with a header notation of '1, 2, etc.'
The finallist is now called finallist.txt instead of shortlist.finallist.txt.
In the log file produced during the StartWithBestScore function SAMPLE_RATIO has been renamed SAMPLE_CV.
The inputDir is not a requirement anymore for the GenerateTargetQCList function.
In the GenerateTargetQC function the searchpattern is changed to prevent the step trying to process files other than normalized coverage files.
The used commandline options and values are now printed to the stdout by default.
Literature
CoNVaDING: Single Exon Variation Detection in Targeted NGS Data
Tutorial
A video tutorial can be found here.
Contact
Mailto:
Lennart Johansson l.johansson@umcg.nl
Freerk van Dijk f.van.dijk02@umcg.nl
Last updated