DeviCNV: detection and visualization of exon-level copy number variants in targeted next-generation sequencing data

Background Targeted next-generation sequencing (NGS) is increasingly being adopted in clinical laboratories for genomic diagnostic tests. Results We developed a new computational method, DeviCNV, intended for the detection of exon-level copy number variants (CNVs) in targeted NGS data. DeviCNV builds linear regression models with bootstrapping for every probe to capture the relationship between read depth of an individual probe and the median of read depth values of all probes in the sample. From the regression models, it estimates the read depth ratio of the observed and predicted read depth with confidence interval for each probe which is applied to a circular binary segmentation (CBS) algorithm to obtain CNV candidates. Then, it assigns confidence scores to those candidates based on the reliability and strength of the CNV signals inferred from the read depth ratios of the probes within them. Finally, it also provides gene-centric plots with confidence levels of CNV candidates for visual inspection. We applied DeviCNV to targeted NGS data generated for newborn screening and demonstrated its ability to detect novel pathogenic CNVs from clinical samples. Conclusions We propose a new pragmatic method for detecting CNVs in targeted NGS data with an intuitive visualization and a systematic method to assign confidence scores for candidate CNVs. Since DeviCNV was developed for use in clinical diagnosis, sensitivity is increased by the detection of exon-level CNVs. Electronic supplementary material The online version of this article (10.1186/s12859-018-2409-6) contains supplementary material, which is available to authorized users.


Note S1. Performance comparison based on the mean target depth for a sample
To identify the minimum mean target depth in samples for detecting real CNVs, we ran DeviCNV with down-sampled samples. We selected one sample in a batch and made 10 input sets by down-sampling it by uisng samtools [1]. Total 40 input sets were created by repeating this process for four samples. We ran DeviCNV using these 40 sets as inputs. We examined whether the known CNVs were found in these down-sampled samples for each input set, and the number of detected candidates was measured.
In NA06804, known CNVs were not found when the mean target depth dropped to below 100X. Similarly, mean target depth dropped to below 60X for NA22208 sample, so no known CNV was detected. In the case of NA01741, the known CNV was well-detected regardless of how low the mean target depth was for the sample. These samples assumed that they have good uniformity and strong CNV signal. Likewise, the lower the mean target depth, the more the number of CNV candidates that were selected. Lastly, in GM14603, which had low uniformity, the number of CNV candidates did not change according to its mean target depth. As a result, this sample showed too many CNV candidates. In the case of IMD_PCR GM14603, similarly to IMD_HYB samples, the read depth of each pool dropped below 100, so no known CNV was detected.  PCR, polymerase chain reaction based capture approach; EX, exon; DEL, deletion; DUP, duplication.
a Indicates whether each method found a known CNV. "O" means all CNVs were found, "X" means they were not found at all, and "A" means they were found only in some exons; b indicates the number of CNV candidates that received the highest score of 5 found in the corresponding sample.
To evaluate the performance of DeviCNV, we performed qPCR on the subset of CNV candidates with confidence score of 5 from the IMD_HYB dataset. The subset was selected from 11 cell lines with the number of CNV candidates of score 5 less than 10, which resulted in a total of 40 CNV candidates (27 duplications and 13 deletions). We randomly selected 25 of the 497 CNV candidates with confidence score of 4 from the above 11 cell lines. We selected 25 CNV candidates in 10 samples.
The a Indicates whether each method found a known CNV. "O" means CNVs were found, and "X" means they were not found at all.

Inclusion criteria
(1) Patients who have been commissioned by primary/secondary hospitals or Seoul National University Children's Hospital newborn room to department of Pediatrics in SNU Children's Hospital for additional confirmation after abnormal findings in Inherited metabolic screening, which is performed on all newborns in Korea.
(2) Patients currently in treatment, diagnosed from secondary blood/urine metabolism testing without identified genetic cause, after abnormal findings in initial IMD screening in Korea.

Exclusion criteria
(1) Low birth weight less than 2.0 kg (2) Premature infants less than 35 weeks

Note S5. Visual inspection process to find pathogenic CNVs in patients
DeviCNV was used to find pathogenic CNVs in patients with inherited metabolic disorders through the IMD panel.
First, a CNV candidate list table and plots (one whole-genome view plot and single-gene view plots) were generated by DeviCNV, and SNV/INDEL analysis results were also obtained using SNV/INDEL caller with patient NGS data. We also selected genes related to patient phenotype by analyzing clinical information of patients. By combining these three sets of information, our clinician selected pathogenic CNV through visual inspection. The order of inspecting CNV candidates was as follows: (1) within phenotype-related gene and with high score, (2) within phenotype-related gene and with low score, and (3) others. To perform the visual inspection, the clinician first looked at the overall CNV calling quality through the whole-genome view plot of the sample. The clinician checked that the mean target depth of the sample was not too low, how many data points were classified as lowquality type, and that there were not too many duplications or deletions. Next, the clinician looked at the gene-centric view plot with the CNV candidate. The clinician identified patterns of data points belonging to the CNV candidate, and confirmed that the expected read depth ratio was sufficiently low or increased.
Through this process, the clinician selected pathogenic CNV candidates and confirmed them by qPCR.

Note S6. Performance comparison based on the number of input samples
CNV, copy number variation; IMD, inherited metabolism disorder; HYB, hybridization based capture approach; EX, exon; DEL, deletion; DUP, duplication; NA, not available.
a Indicates whether each method found a known CNV. "O" means all CNVs were found, and "X" means they were not found at all; b indicates the number of CNV candidates that received the highest score of 5 found in the corresponding sample.

Note S7. Performance comparison based on the configuration of the sample set used as an input
We recommend that the input should consist of a sample set from the same batch. We compared CNV detection performance when using an input set obtained by mixing samples from different batches with when using an input set with samples from the same batch.
We used three batch sets of data from the IMD_HYB: original_batch1, original_batch2 and original_batch3. We randomly mixed 66 samples to construct three random batch sets: random_batch1, random_batch2 and random_batch3. We also combined all the samples into a single set without regard to their batch: "merged all" batch.
Each set was analyzed using DeviCNV_HYB. Of 66 samples, seven samples in the original batch, 10 samples in the random batch and 10 samples in the combined batch were filtered out, due to a low correlation coefficient. For comparisons under the same conditions, the mean and median of the raw and highest confidence CNV candidates were measured, except for the samples that were filtered out in any of the three experiments. In experiments using the original batch, a small number of CNVs was detected. In contrast, the largest number of CNVs was detected from experiments using the random batch and "merged all" batch.
For the "merged all" batch experiment, the number of input samples was three times that of the random batch experiment. Since there were more data to use for modeling, we expected that performance should improve and detect fewer CNV candidates. Alternatively, in both experiments, we expected that CNV candidates should be found in similar numbers since the proportion of samples from the same batch was the same in the input sample set configuration. However, we observed less CNV candidates in the random batch experiment than in the "merged all" batch experiment.
To analyze this, we plotted the median read depth of the sample as X value, and the number of raw CNV candidates found in the "merged all" batch minus the number of raw CNV candidates found in the random batch as the Y value. As a result, we found that the samples with low median read depth (<200X) weirdly got more CNV candidates in the "merged all" batch. Therefore, users need to look carefully at the results when calling the CNV of a sample with depth less than 200. a Indicates whether each method found a known CNV. "O" means all CNVs were found, and "X" means they were not found at all; b indicates the number of CNV candidates that received the highest score of 5 found in the corresponding sample.

Note S8. Differences in the number of data points for each exon based on input intervals
Usually, one more probes tiled in one exon for targeted NGS panel design. When we analyze NGS data, we considered which unit is better to detect CNV signal: individual probes or individual exons (or merged probes).
DeviCNV does not merge overlapping probe (or amplicon) intervals, which would cause reduction in probe-specific information; instead, it analyzes probes separately. Therefore, a probe-level read depth is estimated and a probe-level read depth ratio is predicted, and DeviCNV can obtain one more signals for one exon.
In contrast, if a merged probe interval is used for a single exon or a target interval (usually an exon interval) for one exon, usually only one signal can be collected for each exon. In the case of fewer data points per exon, it is not likely to detect single exon CNVs; and even if some CNVs are detected, they would be less reliable.
To confirm differences of detecting performances between using individual probes and merged probes (=target intervals), we ran DeviCNV to analyze IMD_HYB and IMD_PCR samples with probes and targets(exon) as interval inputs.
In the experiment with probes as the interval input, we used default parameter setting for scoring system. In the experiment with targets as the interval input, we extracted duplication/deletion regions covered by one or more consecutive probe-level significant duplications/deletions at the step of adding duplication or deletion region covered by consecutive strong probe-level CNV signals, and changed the threshold of "ProbeCntInRegion" filter to one or more and removed the threshold of "STDOfReadDepthRatios" filter to detect single-level exon CNVs. With total 14 cell-line known CNVs, we compared the performances of two experiments. In the PCR-based method using multiple pools, it was more useful to use amplicon information rather than target information as it detects fewer CNV candidates. This is thought to be due to low information loss. For the hybridization-based method, if probe information could not be obtained, the target information was used as the input interval.
However, the probe interval input method with less information loss showed better performance when looking for small size CNVs (GAA EX18 DEL). a Indicates whether each method found a known CNV. "O" means all CNVs were found, and "X" means they were not found at all; b indicates the number of CNV candidates found in the corresponding sample. The number of 5-score CNV candidates that received the highest score for the probe interval input method and the number of 4-score CNV candidates that received the highest score for the target interval input method is indicated.

Note S9. Performance comparison based on MQV thresholds
DeviCNV uses only reads that exceed the user-input mapping quality value (MQV) threshold. We compared cell-line CNV deletion results from DeviCNV using the following MQV thresholds: MQV≥0 and MQV≥20.
The statistics of raw CNV candidates, before scoring with CNV candidates by the score system, are shown in the following Average, an average of CNV candidates with the highest score of 5 of all samples except low-quality samples; Median, a median of CNV candidates with the highest score of 5 of all samples except low-quality samples; DEL, deletion; DUP, duplication; TOTAL, the total number of deletions and duplications.
a Indicates how many input samples fail in the low-quality sample filtering step.
Comparison of the two experiments revealed that the number of the highest confidence CNV candidates and the performance of the known cell-line CNV detection were similar.
The MQV threshold may affect the read depth calculation in authentic gene regions for genes with corresponding pseudogenes, since the reads in these regions had low MQV values.
In case of GAA single-exon deletion of GM14603 of IMD_HYB, the coverage of the corresponding sample was not good. Therefore, we confirmed that this CNV was only detected in MQV0 experiment, using as many reads as possible.
For the 1-copy deletion of CYP21A2 of NA12217 samples of IMD_HYB/IMD_PCR, it was too difficult to detect due to the effect of pseudogene. We confirmed that the CNV was detected in the IMD_HYB MQV≥20 experiment using only reliable reads.
For the 1-copy deletion of OTC of GM24007 sample of IMD_PCR, the quality of this sample was bad for CNV calling, due to low coverage and filter-out from low-quality sample filter. Therefore, OTC deletion was not detected in all experiments. a Indicates whether each method found a known CNV. "O" means all CNVs were found, and "X" means they were not found at all; b indicates the number of CNV candidates that received the highest score of 5 found in the corresponding sample.

Note S10. Low-quality sample filter by using sample-to-sample correlation
To filter out low-quality samples, we assumed that if the number of segments resulting from running a CBS (circular binary segmentation) of a sample is too many than other samples, the sample is considered as low-quality, which has high probability to detect abnormally large CNVs. Since large amount of CNV candidates from these low-quality samples are unreliable, they should be excluded when creating a linear regression model.
DeviCNV finds CNV candidates by comparing the read depth ratios of the samples in the same batch. Therefore, for samples with low sample-to-sample correlation values, the number of CNVs candidates must be high. Therefore, these samples must be filtered out in the generating model step.
To determine the criteria for the low-quality sample filter, we analyzed the histogram of CBS segments from the all samples of IMD_HYB and IMD_PCR. To count CBS segments for each sample, only the number of segments with read depth ratio values of more than 1.3 or less than 0.7 (possibly a candidate for duplication/deletion) was selected.
Samples with 250 or more segments in IMD_HYB and samples with 50 or more in IMD_PCR were defined as low-quality samples.
Then, we compared the number of segments of the samples, as well as the sample-tosample correlation with other samples. In a batch of n samples, we compared one sample to the other samples (n-1) to obtain (n-1) correlation values for the 100 th , 75 th , 50 th , 25 th percentiles and the minimum, respectively. Each of these five sample-to-sample correlation values were compared with the number of segments in the sample, and the correlations between them were calculated. Finally, we selected the one with the highest correlation. In the case of IMD_PCR analysis, the association between the number of segments and the sample-to-sample correlation value was low. For the IMD_HYB analysis, the 75 th percentile value of the sample-to-sample correlation was -0.9363, which was most closely related to the number of segments.
Finally, the 75th percentile sample-to-sample correlation value was considered to filter out the sample with a value less than 0.7.

Note S11. Performance comparison based on duplication and deletion thresholds for read depth ratios
When DeviCNV calculates p-values of probe-level CNVs and select raw CNV candidates, it needs duplication and deletion thresholds. We performed experiments with five duplication/deletion thresholds using IMD_HYB and IMD_PCR.
The following table shows statistics for raw CNV candidates before scoring. The stricter the duplication/deletion thresholds, the smaller the number of raw CNV candidates detected. The following table shows statistics for CNV candidates with confidence score of 5, which is the highest score from the scoring system of DeviCNV. DeviCNV extracts reliable CNV candidates through a scoring system. If the duplication/deletion thresholds are too low, it occurs a high probability that the pattern of data points in one CNV candidate is not stable, or the average read depth ratio value is not low(or high) enough. Therefore, it is understandable that 5-score CNV candidates are extracted less in the "TH.dup: 1. a Indicates whether each method found a known CNV. "O" means all CNVs were found, "X" means they were not found at all, and "A" means they were found only in some exons; b indicates the number of CNV candidates that received the highest score of 5 found in the corresponding sample.
Note S12. Unique scoring system for selecting high-confidence CNV candidates (1) ProbeCntInRegion: how many signals support the CNV candidate We count the number of probe-level read depth ratios for the CNV candidate. The larger this value is, the higher the reliability of the CNV candidate, because this indicates that many signals support that candidate. If the user only wants to get CNV candidates supported by more multiple probes, user can increase this value.
(2) AverageOfReadDepthRatios: how strong the signal supporting the CNV candidate This value represents how far an average of the median read depth ratios supporting the CNV candidate is from the normal state (read depth ratio = 1). In the gene-centric plot, the red rectangle indicates the median read depth ratio, and the average of these values is AverageOfReadDepthRatios. The smaller this value is, the higher the reliability of detecting a deletion candidate. In contrast, the larger this value is, the higher the reliability of detecting a duplication candidate.
(  genes, covering 150,802 bps. This is the only version including a PCR-based approach.

Note S14. Generating targeted NGS data
Targeted NGS data from cell lines and clinical samples were generated for validation and identification of known cell-line copy number variants (CNVs) and pathogenic CNVs.
To briefly describe the hybridization-based capture approach, sample DNA was sheared to around 150 bp, and barcoded adapters were ligated to allow the DNA fragments to be captured and sequenced. The fragments belonging to the target region were captured, purified, and sequenced using the 100 bp paired-end setting on a HiSeq instrument.
For the PCR-based capture approach, locus specific primers with barcoded adapters were designed to attach to a targeted region during PCR. Following PCR amplification, the amplified DNA fragments were pooled, purified, and sequenced on an Ion Torrent Personal Genome Machine (PGM) system or S5 (Thermo Fisher Scientific, Waltham, MA, USA). The data were processed for alignment with the standard Ion Torrent Suite™ Software on the Torrent Server.
Note S15. Failure rate of DeviCNV, VisCap, XHMM, and CODEX The performance of DeviCNV was compared to VisCap, XHMM, and CODEX. The following table shows failure rates after sample QC for each tool.
In the case of VisCap, there is 'run_1' which calls CNV using all samples. After that, 'run_2', which calls CNV again, is repeated while removing samples with low quality. This process is performed until all remaining samples pass their own low-quality sample filter, and the number of runs per batch is different. The following table shows failure rates for 'run_1', 'run_2', and 'run_last', the last run of the batch. For most batches in all panels of VisCap, the failure rates were high when 'run_2'. Therefore, the result of 'run_1' was used for comparison with DeviCNV.