Considerations for clinical read alignment and mutational profiling using next-generation sequencing

Introduction

Since emergence in 2005, next-generation sequencing (NGS) technologies have proven prolific tools in the research setting, permeating a variety of scientific disciplines and demonstrating a range of applications that seems to be limited only by the imagination of the sequencing community. The technology continues to develop at a rapid pace with established instrument manufacturers regularly augmenting their product portfolios and an increasing number of start-up companies promising to disrupt the market. Beyond basic research applications, NGS technologies are now increasingly being applied in the clinical environment, driven partly by their rapid maturation and the arrival to market of smaller, cheaper sequencing platforms.

The potential clinical application of NGS has a broad scope ranging from full human genome profiling¹ to investigation of the microbiome² and includes applications such as biomarker discovery, patient diagnosis, prediction of drug response and patient stratification for clinical trials. Such applications often involve the targeted profiling of genes known to be of clinical relevance. These genes harbor diagnostically relevant variants including single nucleotide polymorphisms (SNPs), and small insertions and deletions (INDELs). Individual genes have previously been interrogated in clinical testing using traditional techniques such as Sanger sequencing however NGS technologies have already begun to supplant the previous tools of choice in these areas, offering increased speed and throughput with reduced running costs.

Despite many successes and increasing uptake, the data generated by NGS analyzers is not perfect, with each platform yielding characteristic errors and biases. Furthermore, NGS technologies produce reads that are much shorter than those traditionally produced by Sanger sequencing methods and this can complicate matters further, especially in genomes containing a large proportion of repetitive elements³. The effect of these problems is most visible in large scale studies such as genome-wide sequencing where a recent study reported a 1 million variant platform-based discrepancy for a single genome⁴. This fact bestows responsibility on both algorithm and software developers and downstream users to develop deep understanding of the various data types and their idiosyncrasies and to apply this appreciation in their analysis and interpretation, in order to correct or compensate for potential errors. Despite forming an area of active research, data interpretation remains an issue and is no doubt a factor in feeding the inertia of many clinical facilities that are reluctant to adopt the new technologies⁵.

Two major computational steps in variant detection from NGS data are read alignment whereby the data are mapped to corresponding locations on a target genome, and mutation calling whereby nucleotides differing from the target genome are assessed and scored on their likelihood of representing a genuine mutation versus an error. While these two stages of analysis may be supplemented with various pre or post processing techniques they represent the most crucial steps and therefore the area of most active software and algorithm development.

Aligners of choice have begun to emerge^6,7 however their strengths are often application specific and different tools are recommended depending on the sequencing platform and individual study goals. Aligners generally fall into the categories of being based on either hash tables or suffix trees⁸. Suffix tree-based aligners such as Heng Li’s BWA⁹ are characterized by their speed and memory-efficiency whilst generally achieving lower sensitivity than hash-based counterparts such as Stampy¹⁰. There is thus a necessary trade-off between speed and sensitivity at the read alignment stage with speed often being prioritized due to the volumes of data produced by NGS technologies and the corresponding time required for analysis. Aligners can be further classified as gapped or ungapped based on their ability to produce successful alignments in the presence of small INDELs. Aligners including BWA and Stampy have been shown to produce alignments with a fair degree of success for reads containing a range of INDEL sizes¹⁰ however such abilities will vary based on a range of complicating factors including size and location of the INDEL. As well as generating continuous controversy¹¹, the BRCA1 gene presents a particularly interesting set of alignment challenges due to a disproportionately high concentration of INDELs greater than eight nucleotides in length. In fact, 3% of known deleterious mutations in the Breast Cancer Information Core (BIC) database¹² fall into this category and represent a significant over-representation when compared to a healthy genome. Furthermore the BRCA1 gene contains multiple areas of high shared identity in the form of tandem repeats, posing another difficulty in achieving accurate read mapping. Challenges like this pose particular difficulties in the clinical setting where errors have the potential to translate to misdiagnosis or mistreatment, directly affecting and endangering the lives of patients. No gold-standard clinical alignment tools yet exist and numerous publicized examples of early translational work appear to base their choice of tool on user-friendliness or availability of a graphical-user-interface rather than assessments of performance. We have investigated the performance of a range of popular alignment tools and assessed their ability to accurately detect known mutations. Several of these tools are already being used as components of diagnostic workflows in the clinical setting. Here we present data generated using simulated reads derived from the human BRCA1 gene. Our findings demonstrate the widely varying abilities of common read alignment tools and their impact on downstream variant calling. Furthermore the results suggest a need for careful and thorough evaluation of all tools being used in a particular analysis pipeline through simulation and analysis of data of known constitution.

Methods

Results/discussion

Sensitivity of detection

The greatest overall sensitivities of detection were achieved by Novoalign and Omixon’s Variant Toolkit in paired-end and single-end mode respectively (Figure 1, Table 3). Stampy also enabled highly sensitive mutation detection with Bfast performing least favorably. Sensitivities were also assessed based on the category of mutation (Table 4).

Figure 1. Overall sensitivity of detection.

Graphical representation of the total percentage of mutations detected across of 67 groups of simulated Illumina reads in single-end and paired-end mode.

Table 3. Overall detection sensitivities.

The total percentage of mutations detected by each aligner across the 67 groups of simulated Illumina reads was recorded for single and paired-end run-modes. Each read group contained 13 INDELs of varying sizes as well as 20 SNPs.

Aligner	Sensitivity PE	Sensitivity SE
BWA	99.23	94.53
Bfast	90.68	87.83
Smalt	98.55	98.51
Stampy	99.46	99.05
Mosaik	97.06	95.07
CLC	98.01	97.92
CLCBeta	96.88	96.74
Novoalign	99.59	99.00
Omixon	99.41	99.14
Bowtie2	98.91	98.19
Nextgene	98.64	98.69

Table 4. Detection sensitivities by mutation type.

Total percentage of SNPs and INDELs detected across the 67 groups of simulated Illumina reads by each aligner in single and paired-end mode.

Aligner	% SNPs found PE	% SNPs found SE	% Insertions found PE	% Insertions found SE	% Deletions found PE	% Deletions found SE
BWA	99.48	99.48	99.38	89.38	98.55	85.48
Bfast	92.69	89.85	83.13	78.75	90.20	88.20
Smalt	99.40	99.40	96.56	96.25	97.64	97.64
Stampy	99.48	99.25	99.38	98.13	99.46	99.09
Mosaik	97.01	96.34	96.25	90.00	97.64	94.92
CLC	99.40	99.40	96.56	96.25	95.46	95.28
CLCBeta	97.61	97.61	95.31	95.00	96.01	95.64
Novoalign	99.70	99.48	100.00	98.13	99.09	98.37
Omixon	99.40	98.96	99.38	99.38	99.46	99.46
Bowtie2	99.50	99.50	97.50	96.60	98.40	96.00
Nextgene	99.00	99.00	96.60	96.60	99.10	99.10

Some aligners such as Smalt, Bowtie 2 and CLC were clearly seen to perform more strongly on detection of SNPs than INDELs in general. Nextgene performed similarly for SNPs and deletions but had lower sensitivity of insertion detection. BWA showed obvious decreases in ability to detect INDELs when shifting from paired-end to single-end mode. In contrast, Novoalign, Omixon and Stampy performed well regardless of run-mode or mutation type. Overall Novoalign was the best performer in paired-end mode while Omixon achieved the highest sensitivities in single-end mode. In a clinical environment, sensitivity will likely represent the most important metric in evaluating alignment software, however other factors may also be of importance, depending on the test in question. Specificity values are not included here as they were non-discriminatory in this context due to low numbers of false positives relative to the high number of true negatives.

Incorrect identification of mutations

Controlling the rate of false positive results is clinically important in avoiding unnecessary treatment, expense and patient anxiety. The number of incorrectly identified mutations varied widely dependent on aligner and run-mode (Figure 2, Table 5). The highest rates were obtained using Bfast, followed by CLC Bio’s Beta aligner, Nextgene, Mosaik, Stampy, and Bowtie 2 while Novoalign and Smalt were the strongest performers followed closely by Omixon and CLC. Number of false positives alone is of limited utility in assessing aligner performance, however positive predictive value (PPV) provides a useful metric which combines counts of both true and false positives in a single value (Table 6). PPV was calculated based on the equation:

$$\mathit{\text{PPV}}=\frac{\mathit{\text{True Positives}}}{\mathit{\text{True Positives}}+\mathit{\text{False Positives}}}$$

Figure 2. Graphical representation of the total number of false positives expressed as a percentage of the true positive mutations detected by each aligner in single and paired-end mode across 67 groups of simulated Illumina reads.

Each read group contained 20 SNPs and 13 INDELs.

Table 5. Raw numbers of false positive and negative mutations detected across 67 groups of simulated Illumina reads by all aligners in single and paired-end mode.

Each read group contained 20 SNPs and 13 INDELs.

Aligner	# False positives PE	# False positives SE	# False negatives PE	# False negatives SE
BWA	23	168	17	121
Bfast	884	759	206	269
Smalt	10	14	32	33
Stampy	263	95	12	21
Mosaik	314	76	65	109
GLG	21	23	44	46
GLGBeta	119	566	69	72
Novoalign	6	10	9	22
Omixon	39	21	13	19
Bowtie2	249	66	24	40
Nextgene	395	117	30	29

Table 6. Positive predictive values (PPVs) calculated for each aligner in single and paired-end mode based upon detection of mutations across 67 groups of simulated Illumina reads each containing 20 SNPs and 13 INDELs.

In this case PPV is the quotient of the total number of true positives divided by the sum of the total number of true positives and true negatives.

Aligner	PPV SE	PPV SE
BWA	98.96	92.56
Bfast	69.40	71.90
Smalt	99.54	99.36
Stampy	89.32	95.84
Mosaik	87.24	96.51
CLC	99.04	98.95
CLCBeta	94.74	79.08
Novoalign	99.73	99.55
Omixon	98.26	99.05
Bowtie2	89.78	97.05
Nextgene	84.67	94.91

No inferences were made about prevalence in calculating the values. Novoalign had the highest PPV in both single and paired-end mode. Smalt, CLC and Omixon’s Variant Toolkit also performed strongly on this metric. Notably Stampy performed relatively poorly in contrast to its high sensitivity.

Effect of INDEL size on detection

Not unexpectedly there appeared to be a general trend of INDEL size affecting most aligners ability to facilitate downstream calling (Figure 3). The size of the effect varied by aligner with BWA and Novoalign showing good detection rates for all but the largest deletions while others such as Bowtie 2 and the CLC aligners were not successful far beyond a 10bp INDEL size. Nextgene showed better sensitivity for insertions than deletions. Stampy and Omixon’s Variant Toolkit were the only two aligners to detect the largest deletions. This highlights a need for those involved with testing and analysis to develop an appreciation of the various mutations that might exist in their target genes and to select their analysis tools appropriately. INDELs in the range represented by the BRCA1 mutation panel have real-world relevance in genetic disorders and strong aligner performance on larger INDELs appears to be an exception rather than a rule.

Figure 3. Effects of INDEL size on success of detection for all aligners in single-end (top panel) and paired-end (bottom panel) mode across 67 groups of simulated Illumina reads, each containing 13 INDELs of varying sizes.

Summary and conclusions

Using a simulated, targeted sequencing scenario with Illumina read data, the work presented here highlights several important considerations regarding aligner choice in studies involving profiling of mutations. Furthermore the data presented goes some way to characterizing the performance of a comprehensive selection of commonly used aligners and should represent a useful resource for anyone focused on similar scientific studies.

Whilst the dataset used in this study was engineered to include a challenging range of mutations and efforts were made to simulate the error profile of Illumina sequencing technology, it is nevertheless a simplified representation of real-world data. Experimental artifacts such as PCR stutter have the potential to present further challenges to alignment algorithms and there is no consideration of such issues here. Furthermore, the test dataset used in this study produced a uniform, high coverage of the target gene and only homozygous variants were simulated. Finally, the use of only a single variant-caller in the study means that some of the errors encountered may not be due to alignment issues. The aim is to follow up the current study by focusing on an expanded gene-set, alternative variant-callers, homo and heterozygous mutations and different sequence formats. Nonetheless, this focused study demonstrates the utility of simulated data in assessing program performance.

With the exception of Bfast, all aligners performed relatively well on the BRCA1 dataset. Clinical applications necessitate the use of the most highly accurate solutions, however. Only Novoalign, Omixon’s Variant Toolkit and Stampy achieved 99% or greater sensitivity in both paired-end and single-end mode. Omixon and Stampy were the only two aligners to detect the longest deletions in the dataset however Stampy’s performance was let down by a false positive rate which would be considered unacceptably high for most applications. While Novoalign did not detect the largest deletions, it was the fastest aligner and the most sensitive in paired-end mode. Nevertheless, assuming the longer run-times are not an issue, Omixon’s superior sensitivity in single-end read mode likely makes it the best option when paired-end protocols are not possible. It should also be mentioned that as a freely available, open-source tool, BWA’s paired-end performance is laudable and goes some way to justifying its widespread use and popularity.

While the tests here produce some clear winners, they also serve to highlight that program performance can vary widely based on the fine-details of a particular run. Even the strongest overall performer can be found lacking in some respects and this means that researchers should be vigilant in their selection of tools for a particular application. In certain instances it may even be necessary to combine two or more approaches to ensure that all relevant aspects of a given dataset are sufficiently characterized and any approach will still require some level of visual inspection and quality control in a clinical setting. The data presented here should facilitate and expedite selection of the correct aligner for a particular task but they do not obviate the requirement for careful consideration nor further testing and analysis on the part of the end-user.