Comparison of somatic mutation calling methods in amplicon and whole exome sequence data
© Xu et al.; licensee BioMed Central Ltd. 2014
Received: 20 December 2013
Accepted: 20 March 2014
Published: 28 March 2014
Skip to main content
© Xu et al.; licensee BioMed Central Ltd. 2014
Received: 20 December 2013
Accepted: 20 March 2014
Published: 28 March 2014
High-throughput sequencing is rapidly becoming common practice in clinical diagnosis and cancer research. Many algorithms have been developed for somatic single nucleotide variant (SNV) detection in matched tumor-normal DNA sequencing. Although numerous studies have compared the performance of various algorithms on exome data, there has not yet been a systematic evaluation using PCR-enriched amplicon data with a range of variant allele fractions. The recently developed gold standard variant set for the reference individual NA12878 by the NIST-led “Genome in a Bottle” Consortium (NIST-GIAB) provides a good resource to evaluate admixtures with various SNV fractions.
Using the NIST-GIAB gold standard, we compared the performance of five popular somatic SNV calling algorithms (GATK UnifiedGenotyper followed by simple subtraction, MuTect, Strelka, SomaticSniper and VarScan2) for matched tumor-normal amplicon and exome sequencing data.
We demonstrated that the five commonly used somatic SNV calling methods are applicable to both targeted amplicon and exome sequencing data. However, the sensitivities of these methods vary based on the allelic fraction of the mutation in the tumor sample. Our analysis can assist researchers in choosing a somatic SNV calling method suitable for their specific needs.
Somatic point mutation calling from matched tumor-normal patient samples is a critical step for cancer genome characterization and clinical genotyping [1, 2]. Over the last few years, next-generation sequencing (NGS) has become a popular strategy for genotyping, enabling more precise mutation detection compared to traditional methods due to its high resolution and high throughput . Nevertheless, detecting somatic mutations is still challenging, especially for low-allelic-fraction variants caused by tumor heterogeneity, copy number alteration, and sample degradation [4–6]. Errors in base calling and read alignment present additional challenges for achieving sensitive and specific somatic variant calls. In recent years, several methods have been developed to enhance somatic mutation calling accuracy [7–11]. In general, these methods belong to two families: (1) independent analysis for tumor and normal datasets from an individual followed by SNV type classification using a statistical significance test or simple subtraction (e.g. ); (2) simultaneous analysis for matched tumor and normal datasets using joint probability-based statistical approaches, e.g. SomaticSniper  and Strelka . In general, agreement among different algorithms is relatively low [13, 14], making selection of candidate SNVs for further validation difficult. This disagreement is likely partially due to different error models and prior assumptions underlying each algorithm.
In light of these challenges, there is a need to systematically evaluate the performance of different SNV calling algorithms to help guide best practices in the research community. Recently, several independent studies compared somatic SNV calling tools [15–17]. For example, Roberts et al.  compared the performance of VarScan, SomaticSniper, JointSNVMix2, and Strelka in detecting SNVs from matched cancer-normal Illumina exome sequencing of a chronic myeloid leukemia patient. They reported a large difference in called variant sets among algorithms. In a separate study, Wang et al.  evaluated 6 tools (EBCall, JointSNVMix, MuTect, SomaticSniper, Strelka, and VarScan2) on tumor-normal pairs using both whole-genome sequencing and exome sequencing data. According to their evaluation, MuTect achieved the highest sensitivity for calling low-allelic-fraction somatic variants. Most of these studies utilized publically available exome data from hybridization capture enrichment, and some used additional validation data from alternative technologies. Complimentarily, Spencer et al.  assessed SAMtools, Genome Analysis Toolkit (GATK), VarScan2, and SPLINTER in detecting low-allelic-fraction variants from synthetic mixtures of targeted hybridization capture sequencing reads. They reported that VarScan2 achieved the highest sensitivity but with a high false positive rate, whereas SPLINTER yielded both high sensitivity and specificity.
Exome capture sequencing is typically limited to ~60x read depth (due to cost), and therefore using exome reads limits the potential to detect low-allelic-fraction SNVs in the sample, especially those at less than 5%. Detecting low-allelic-fraction SNVs in the sample is often important for early diagnosis, prevention of drug resistance, and detection of residual tumors. Targeted panel sequencing has been widely used to increase SNV detection sensitivity by achieving a much higher median read depth (>500x), even on a bench top sequencer [18, 19]. Among the different enrichment methods available for targeted sequencing, PCR amplicon-based enrichment has been widely popular for small- to medium-sized gene panels due to its low input requirement, simple protocol, fast turnaround time, and good performance on FFPE samples. While there are differences between read properties from PCR amplicon sequencing and exome sequencing, many of which could affect SNV calling, there has not yet been a detailed comparison of SNV calling algorithms using reads from PCR amplicon sequencing. In addition to higher read depth, reads from PCR amplicons have more uniformly defined ends, more PCR duplication, and may contain more amplification bias. Because many variant calling algorithms were developed using whole genome or exome data, the question remains of whether performance differs when these algorithms are applied to amplicon sequencing data.
Another challenge in comparing variant calling methods is to identify a set of true positive and true negative variants in the sample, ideally at various allele fractions. This task is not straightforward, due to the lack of sufficient available validation data. One standard assessment method is to check the overlap between variant sets as a way to assess the concordance and discrepancies among different calling methods. Alternatively, some researchers have randomly chosen a few candidates for Sanger sequencing or used SNP arrays for validation. However, Sanger sequencing involves high cost and labor, making the generation of large-scale validation data difficult. In addition, some of the low-allelic-fraction SNVs are very difficult to validate by Sanger sequencing or SNP array due to the detection sensitivity of these technologies (typically insufficient for variants below 20% fraction). To overcome these limitations, the NIST-led “Genome in a Bottle” Consortium (NIST-GIAB)  developed a community resource of high-confidence variants for the reference individual NA12878. From the July 17, 2013 data release, the NIST-GIAB standard is an integration of 12 datasets generated using 5 different sequencing platforms, with the reads aligned using 7 different aligners, and variants called using 3 different variant callers. The availability of the NA12878 reference sample (from Coriell Institute) along with a high-confidence variant set facilitates empirical evaluation of all steps in NGS analysis: enrichment, library preparation, sequencing, alignment, and variant calling. The NIST-GIAB NA12878 variant standard is the basis of emerging tools such as GCAT (http://www.bioplanet.com/gcat), which provides a community platform for comparing variant calling approaches.
In order to evaluate the performance of popular SNV callers on targeted amplicon sequencing reads, especially for low-allelic-fraction SNVs, we mixed reference DNA from NA12878 with DNA from another 1,000 Genomes Project sample (NA19129) at various ratios, thus creating virtual tumor-normal samples with NA12878 “tumor” variant frequencies ranging from 4% to 100%. We then subjected the dilution series samples to multiplex PCR-based enrichment and Illumina MiSeq sequencing. Here we present a comparison of 5 somatic SNV callers using targeted amplicon sequencing from these virtual tumor-normal samples against the NIST-GIAB gold standard variant set. To explore possible performance differences of the SNV callers on amplicon versus exome sequencing data, we also did a similar in silico dilution series analysis using publicly available exome data from NA12878. Our results demonstrate the strengths and weaknesses of each method in calling somatic SNVs at diverse levels of allelic fraction, and also the value of the newly established NIST-GIAB gold standard in facilitating such comparisons.
List of somatic SNV calling methods and versions
Because most matched tumor-normal variant callers were developed using exome sequencing data, we next assessed whether their respective performance differs between high-read-depth PCR-based amplicon sequencing and low-read-depth hybridization capture-based exome sequencing. The two approaches have different read depths and amplification bias, which could affect variant calling accuracy. Therefore, we compared the abilities of the aforementioned five methods to accurately call somatic SNVs in exome sequencing data. For exome reads, we chose a different Yoruba individual (NA18489) than we used for amplicon sequencing (NA19129), because exome sequencing reads for NA18489 were generated using the same capture procedure used for NA12878. The exome sequencing data of NA12878 and NA18489 are publicly available, as part of the 1,000 Genomes Project (Additional file 1: Table S2). We downsampled the preprocessed BAM files and, in silico, mixed the resulting subsets of reads to desired “tumor” sample purities. The resulting in silico mixture series is diluted in the same way as the amplicon sequencing experimental dilution series: at NA12878 unique heterozygous (homozygous alternative) variant positions over the region of interest, the median variant allele fraction is 4% (8%), 8% (16%), 18% (36%) and 50% (100%), respectively (Figure 1B). Then we applied the five variant calling methods on the in silico matched “tumor-normal” exome sequencing pairs over the region of interest (Additional file 1: Table S1). Additionally, we checked performance results from the CCP region of the exome sequencing data and verified that there are no special advantages/disadvantages to the CCP region relative to the exome capture region (data not shown).
Overall, the sensitivity ranking of each method is comparable with what was achieved in targeted amplicon sequencing data (Figure 2). For the exome data, MuTect and Strelka detected more somatic mutations than the other methods. MuTect gave the best sensitivity among lower concentrations and Strelka gave the second-best sensitivity. Results from the MuTect publication  show that for variants at 5% abundance, MuTect’s sensitivity was between 0.16 and 0.53 for sequencing depth between 30x and 60x, whereas Strelka, used with standard settings, detected between 0.16 and 0.40 in the same sequencing depth range. Therefore, our present observations are consistent with this previous study. SomaticSniper, VarScan2, and NaiveSubtract still yielded low sensitivity in detecting lower-allelic-fraction somatic SNVs. In addition, consistent with what was observed in the targeted sequencing data, all methods achieved high sensitivity in the high-allelic-fraction tumor samples.
We next assessed specificity using only the C50 sample in order not to introduce potential artifacts arising from the in silico mixing. MuTect and Strelka achieved the highest specificity (Figure 3). All methods except NaiveSubtract showed a lower false positive rate compared to amplicon sequencing. This is could be partially due to the fact that all the methods were primarily designed on genome-wide sequencing data or exome data and therefore better control false positives in that data type. Furthermore, since false positive rate is positively correlated with sequencing read depth , the greater number of reads generated by amplicon sequencing technology may yield more false positives.
As the first step in analyzing cancer sequencing data, detecting variants with high sensitivity and specificity is of great importance. Our results highlight the differences in performance between five algorithms when used with targeted amplicon versus exome sequencing data. When analyzing different types of data, use of different algorithms may be appropriate. For example, methods using a hypothesis of Poisson/binomial distribution of alternate allele read depth would deviate from actual exome sequencing data, which display a notable overdispersion in the distribution of alternate allele read depth . In terms of sensitivity, we observed that the five algorithms perform better in calling low-allelic-fraction variants from amplicon sequencing data than from exome data, probably due to the higher sequencing depth in the former. In fact, when the amplicon sequencing data were in silico downsampled to the same sequencing depth as the exome data, the sensitivity from all algorithms dropped and declined more significantly in C04 and C08, while the relative performance ranking remained the same (data not shown). On the other hand, since sequencing read depth is among the factors affecting false positive rate, most algorithms generate more false positives in amplicon sequencing data. However, we also demonstrate through Strelka and MuTect that further tuning the “joint probability” score away from the default setting can reduce the false positives in amplicon sequencing to levels similar to that for exome sequencing, without substantially compromising sensitivity. In summary, the relative performance ranking among five variant callers examined in hybridization capture-based exome sequencing largely holds true for PCR-based amplicon sequencing. Amplicon sequencing can achieve much higher variant detection sensitivity for low-allelic-fraction SNVs. At the same time, default variant calling settings should be further tuned to dampen false positive rates when analyzing amplicon sequencing data.
Regarding individual algorithms, MuTect and Strelka show similarly high sensitivity in detecting low-allelic-fraction SNVs in amplicon sequencing reads, with good specificity. For high-allelic-fraction somatic SNVs, VarScan2 and SomaticSniper perform slightly better than MuTect and Strelka in terms of sensitivity in both amplicon sequencing and exome sequencing data.
NaiveSubtract is very straightforward in concept. Using NaiveSubtract, variants detected in the tumor but not the normal sample would be considered putative somatic mutations. Early studies applied such a subtraction method in detecting somatic mutations and demonstrated success to some extent . One drawback of NaiveSubtract is that it does not take into account the complexity of germline false positives due to PCR-based enrichment artifacts or sequencing errors, or false negatives due to under-sampling during enrichment. Also, NaiveSubtract is not tailored to handle detection of allele sites with low fraction. In other words, the sensitivity is limited by the actual variant caller used, which, in this case, is GATK UnifiedGenotyper.
The performance of VarScan2 is closely related to the minimum allele fraction threshold, because VarScan2 will suppress any mutations below this threshold. We could improve VarScan2 sensitivity after tuning the minimum allele fraction cut-off from the default 0.20 to a lower value, but at the cost of significantly compromised specificity (data not shown). In one recent study , VarScan2 achieved a sensitivity of 0.97 for detecting variants of 5% fraction in targeted sequencing with average coverage of >1000×. However, the false positive rate in this study was 59 per Mb after the poor quality, low coverage, and germline sites were masked, which is about 3-fold higher than our results. Furthermore, prior knowledge about tumor purity is unavailable in real-world applications, so users often use default settings to run VarScan2. Using the default VarScan2 settings, we were unable to call somatic SNVs of low allelic fraction.
Although we set the parameters to maximize sensitivity, Somatic Sniper’s performance still falls behind other algorithms in calling somatic SNVs of low allelic fraction. However, SomaticSniper did achieve the highest sensitivity in the 100% ‘pure tumor’-normal match. The low sensitivity of VarScan2 and SomaticSniper in calling low-allelic-fraction somatic SNVs is consistent with the results in , which showed sensitivity less than 0.10 for somatic SNVs of fraction lower than 0.20. In contrast, MuTect and Strelka produce high-confidence somatic SNV candidates for any level of allele fraction: MuTect and Strelka can robustly detect somatic SNVs at low allelic fraction, with an average sensitivity of 0.89 and 0.86 for sites above 0.05 allelic fraction, respectively.
The merits of using artificial matched tumor-normal sample mixture series data sets created for calling NA12878 “somatic” SNVs are at least two-fold: (1) the NIST-GIAB gold-standard is a well-defined, high-confidence, large set of SNVs for benchmark evaluation, and (2) this setup allowed us to assess variant calling for somatic SNVs of a wide range of allele fractions from real data. The advantages of the NIST-GIAB gold standard lie in the fact that more sites could be assessed as compared to SNP arrays, which are normally limited to common variants. However, one limitation of our analysis is that a high-confidence variant list for the background germline variants is not available. Therefore, using NA12878 for an artificial “somatic” sample is likely more accurate in evaluating sensitivity than in evaluating specificity. Another limitation is that the NA12878 DNA purchased from Coriell is not exactly the same batch used for the NIST-GIAB variant calls. In addition, one constraint of our artificial matched normal samples is the unmatched error distribution between the “normal” and “tumor” samples. So the virtual tumor-normal matches would potentially reduce the performance of methods relative to performance on real tumor samples that have error distributions more similar to the matched normal sample.
In principle, characterizing how performance of different methods depends on factors such as allele fraction is necessary for designing experiments. The complexity and systematic noise from sequencing from tumor-normal matched samples affect each algorithm, but does so in different ways. For systematic comparison, here we used the same aligned read file for all the methods, whereas impact of aligners on the methods would further distinguish their performance. Deeper understanding of the sample properties and the ultimate goals of a cancer research project will make a difference in the choice of algorithms. This work can be extended to address questions regarding detecting variants from small amount of target genome mixed in the control samples. Our work provides insights for the community to assist in choosing sensitive and specific methods for somatic SNV calling, given a specific set of data type and experimental conditions.
In this work, we demonstrated that the five somatic SNV calling methods are applicable to both matched tumor-normal targeted amplicon and exome sequencing data. Sensitivities vary with allelic fraction of the mutation in the tumor sample. Our results show MuTect and Strelka achieved the highest sensitivities. Our analysis can assist researchers in choosing a somatic SNV calling method for their specific needs.
DNA samples of NA12878 and NA19129 were purchased from Coriell Institute. Sample mixtures were created based on the actual amplifiable DNA in each sample, resulting in 0%, 8%, 16%, 36%, and 100% of NA12878 sample mixed in the NA19129 sample, respectively. We treated the mixed samples at 8%, 16%, 36%, and 100% as the virtual tumor samples and the 0% as the virtual normal sample. Therefore, we produced artificially matched tumor-normal samples to call somatic SNVs. QIAGEN’s GeneRead DNAseq Comprehensive Cancer Gene Panel (Version 1) was used to amplify the target region of interest (124 genes, 800 Kb). Targeted enrichment and library construction were done following manufacturer’s user manual (QIAGEN). Briefly, 160 ng of each virtual tumor sample or virtual normal sample were PCR amplified and purified. For each sample, one quarter of the PCR product was used for constructing a barcoded Illumina DNA library. Libraries were quantified using QIAGEN’s GeneRead DNAseq Library Quant System, and were mixed in equal amount. Illumina MiSeq sequencing was performed following manufacturer’s user manual (Illumina) to generate FASTQ files, resulting in median coverage depth of 500–700x for each sample. Experiments were performed in triplicate. The sequencing reads were aligned to the reference genome using BWA , then preprocessed according to Broad best practice guidance (indel realignment, base quality score recalibration, and base alignment quality scoring). Finally, the primer sequences were trimmed away from the pre-processed reads before variant calling.
Whole exome sequencing data of NA12878 and NA18489 were obtained from the 1,000 Genomes Project . Raw reads were mapped to the hg19 human reference genome sequence (GRCh37) by the 1,000 Genomes Project. The publicly available read alignment files in BAM format were generated with BWA mapping, GATK local realignment around known indels, GATK base quality score recalibration, and Picard MarkDuplicates. Each BAM file was randomly downsampled using Picard command-line tools with various percentages to achieve the following mixture concentrations: 8%, 16%, 36%, and 50%. One NA12878 downsampled result and one NA18489 downsampled result were then merged using SAMtools . This step was repeated for different mixture percentages to create an in silico dilution series. Consequently, the primary NA18489 sequencing results represent the virtual normal sample and the NA12878 dilution series were treated as the matched virtual tumor samples.
For targeted amplicon sequencing data, the region of interest (ROI) includes 661,684 bases, resulting from the intersection between QIAGEN CCP panel regions and NIST-GIAB high-confidence callable regions (Version 2.15, excluding genomic sites with conflicting genotype evidence, repetitive regions, decoy sequences, or reported structural variants, as well as regions/variant locations with evidence of bias (systematic sequencing errors, local alignment bias, etc.)), plus masking 213 bp of NA19129 variant sites annotated by 1,000 Genomes Project (ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/release/20110521/). For exome sequencing data, the region of interest include 33,886,321 bases, resulting from the intersection between the exome capture regions with at least one read coverage in both samples and the NIST-GIAB high-confidence sites, plus masking 12,974 bp of NA18489 variant sites annotated by 1,000 Genomes Project (ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/release/20110521/). For these regions of interest, the NIST-GIAB gold standard contains a 181 and 9,868 ‘somatic’ SNVs for our targeted-sequencing and exome-sequencing evaluation, respectively (for details see Additional file 1: Table S1).
NaiveSubtract — SNVs were called separately from virtual tumor and normal samples using GATK UnifiedGenotyper . For exome sequencing data, reads were already mapped, locally realigned and recalibrated by the 1,000 Genomes Project. So SNVs were directly called on the BAM files using GATK Unified Genotyper. Then, SNVs detected in the virtual normal sample were removed from the list of SNVs detected in the virtual tumor sample, leaving the “somatic” SNVs.
MuTect — MuTect is a method developed for detecting the most likely somatic point mutations in NGS data using a Bayesian classifier approach. The method includes pre-processing aligned reads separately in tumor and normal samples and post-processing resulting variants by applying an additional set of filters. We ran MuTect under the High-Confidence mode with its default parameter settings. We disabled the “Clustered position” filter and the “dbSNP filter” for the amplicon sequencing reads, and we disabled the “dbSNP filter” for the exome sequencing.
SomaticSniper — SomaticSniper calculates the Bayesian posterior probability of each possible joint genotype across the normal and cancer samples. We tuned the software’s parameters to increase sensitivity and then filtered raw results using a Somatic Score cut-off of 20 to improve specificity.
Strelka — Strelka reports the most likely genotype for tumor and normal samples based on a Bayesian probability model. Post-calling filters built into the software are based on factors such as read depth, mismatches, and overlap with indels. We skipped depth filtration for exome and amplicon sequencing data as recommended by the Strelka authors. For the amplicon sequencing reads, we set the minimum MAPQ score at 17 for consistency with the defaults in GATK UnifiedGenotyper. We used variants passing Strelka post-calling filters for analysis.
VarScan2 — VarScan2 performs analyses independently on pileup files from the tumor and normal samples to heuristically call a genotype at positions achieving certain thresholds of coverage and quality. Then, sites of the genotypes not matched in tumor and normal samples are classified into somatic, germline, or ambiguous groups using Fisher’s exact test. We generated the pileup files using SAMtools mpileup command.
The compatibility of the output VCF files between different methods as well as the NIST-GIAB gold standard was examined using bcbio.variation tools and manual inspection. The reported SNP call representations between files are comparable to each other.
The command line scripts of variant calling algorithms we used in the paper are included in the additional files.
We acknowledge QIAGEN’s bioinformatics team member Sandeep Namburi for technical help and technical writer Allison Bierly for manuscript review.
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.