Detailed simulation of cancer exome sequencing data reveals differences and common limitations of variant callers

Abstract

Background: Next-generation sequencing of matched tumor and normal biopsy pairs has become a technology of paramount importance for precision cancer treatment. Sequencing costs have dropped tremendously, allowing the sequencing of the whole exome of tumors for just a fraction of the total treatment costs. However, clinicians and scientists cannot take full advantage of the generated data because the accuracy of analysis pipelines is limited. This particularly concerns the reliable identification of subclonal mutations in a cancer tissue sample with very low frequencies, which may be clinically relevant.

Results: Using simulations based on kidney tumor data, we compared the performance of nine state-of-the-art variant callers, namely deepSNV, GATK HaplotypeCaller, GATK UnifiedGenotyper, JointSNVMix2, MuTect, SAMtools, SiNVICT, SomaticSniper, and VarScan2. The comparison was done as a function of variant allele frequencies and coverage. Our analysis revealed that deepSNV and JointSNVMix2 perform very well, especially in the low-frequency range. We attributed false positive and false negative calls of the nine tools to specific error sources and assigned them to processing steps of the pipeline. All of these errors can be expected to occur in real data sets. We found that modifying certain steps of the pipeline or parameters of the tools can lead to substantial improvements in performance. Furthermore, a novel integration strategy that combines the ranks of the variants yielded the best performance. More precisely, the rank-combination of deepSNV, JointSNVMix2, MuTect, SiNVICT and VarScan2 reached a sensitivity of 78% when fixing the precision at 90%, and outperformed all individual tools, where the maximum sensitivity was 71% with the same precision.

Conclusions: The choice of well-performing tools for alignment and variant calling is crucial for the correct interpretation of exome sequencing data obtained from mixed samples, and common pipelines are suboptimal. We were able to relate observed substantial differences in performance to the underlying statistical models of the tools, and to pinpoint the error sources of false positive and false negative calls. These findings might inspire new software developments that improve exome sequencing pipelines and further the field of precision cancer treatment.

Keywords: Cancer genomics; Exome sequencing; SNV; Variant caller integration; Variant calling.

Fig. 1

Workflow of the comparison of the nine variant callers. A heterogeneous cancer sample is simulated based on a real renal cell carcinoma sample (steps 1, 2 and 3). Two different alignment settings and eight different coverage and normal contamination levels are employed (steps 4 and 5). The variant callers deepSNV, GATK HP, GATK UG, JointSNVMix2, MuTect, SAMtools, SiNVICT, somaticSniper and VarScan2 are run on all bam files (step 6). The performance of the different tools is evaluated and analyzed in detail (step 7). The tools are also assessed when using various pipeline or parameter modifications as described in Section “Pipeline and parameter improvements”. A more detailed description of the pipeline and the evaluation procedure can be found in the Methods Section as well as in Additional file 1: Section B, and Additional file 1: Section C

Fig. 2

Performance comparison of variant callers with default parameters. a Sensitivity of variant callers as a function of the variant allele frequency. To make the predictions comparable we selected the largest set of variants from the top of the list of each caller such that the false discovery rate is still below α. We show plots for α equal to 0.05 (solid lines), and 0.1 (dashed lines). If the tool has a very good precision, the two curves for the two α cutoffs are identical, as it is the case for MuTect and VarScan2. The details on how the sensitivities are displayed can be found in Additional file 1: Section G. b Area under precision recall curve as a function of the coverage. Again the two cutoffs for the false discovery rate α=0.05, and α=0.1 are chosen (see Additional file 1: Section C). The coverages correspond to the five different levels (12, 25, 50, 75, and 100%) displayed in Fig. 1 in step 5

Fig. 3

Categories of variant calling errors depending on the quality of the alignment. The top panel (a and b) shows the error categories for the high confidence false positives (prediction sets with at least 95% precision). The bottom panel (c and d) shows the error categories for the high-frequency false negatives (ground truth allele frequency ≥25%). The left panel (a and c) displays the error sources when running default bowtie2 alignments, and the right panel (b and d) displays the error sources when running more sensitive alignments, which were performed with parameters —very-sensitive -k 20, and then choosing the primary alignment for each read with several alignments (samtools view -F 256), i.e. the “best” option. The definition of the categories can be found in Section “Analysis of error sources”

Fig. 4

The effect of pipeline modifications, parameter changes, and combination strategies. We show the sensitivity for the prediction set with at least 90% precision. a Performance when applying local realignment around indels or the binomial test as a germline filter. b Performance of deepSNV, JointSNVMix2, SAMtools, and VarScan2 with different choices of parameters. Additional file 1: Figure H depicts the performance for all parameters that were assessed. c Performance of rank-combinations and intersections of calls from several tools. From each tool, we took the best version. In particular, deepSNV and MuTect with the binomial test as germline filter, SAMtools with option -C 200, SiNVICT with –qscore-cutoff 60, VarScan2 with the parameter –min-var-freq 0.02, as well as the default runs from GATK HP, GATK UG, JointSNVMix2, and somaticSniper. d Summary barplot displaying the performance of the three best rank-combinations as a comparison to each tool individually. If a tool parameter or pipeline change has been used in the rank-combinations, also the performance of the tool in default mode is shown. The y-axis measures the area under precision-recall curve when allowing a false discovery rate of up to 10% (see Additional file 1: Section C)