Comparison of Short-Read Sequence Aligners Indicates Strengths and Weaknesses for Biologists to Consider

Abstract

Aligning short-read sequences is the foundational step to most genomic and transcriptomic analyses, but not all tools perform equally, and choosing among the growing body of available tools can be daunting. Here, in order to increase awareness in the research community, we discuss the merits of common algorithms and programs in a way that should be approachable to biologists with limited experience in bioinformatics. We will only in passing consider the effects of data cleanup, a precursor analysis to most alignment tools, and no consideration will be given to downstream processing of the aligned fragments. To compare aligners [Bowtie2, Burrows Wheeler Aligner (BWA), HISAT2, MUMmer4, STAR, and TopHat2], an RNA-seq dataset was used containing data from 48 geographically distinct samples of the grapevine powdery mildew fungus Erysiphe necator. Based on alignment rate and gene coverage, all aligners performed well with the exception of TopHat2, which HISAT2 superseded. BWA perhaps had the best performance in these metrics, except for longer transcripts (>500 bp) for which HISAT2 and STAR performed well. HISAT2 was ~3-fold faster than the next fastest aligner in runtime, which we consider a secondary factor in most alignments. At the end, this direct comparison of commonly used aligners illustrates key considerations when choosing which tool to use for the specific sequencing data and objectives. No single tool meets all needs for every user, and there are many quality aligners available.

Keywords: accuracy; alignment; comparison; runtime; short-read sequencing.

Figure 1

Creating an FM-index of the word “knickknack.” The first step is to generate all rotations of the reference sequence and sort them lexicographically. The last column is stored as the Burrows-Wheeler transform (BWT) in (A) and the corresponding suffix array in (B). A rank table is created from the BWT, which lists the occurrence and order of each unique character shown in (C). (D) shows the lookup table, which lists the index of the first occurrence of each character from the first column of the sorted matrix.

Figure 2

Creating a suffix array and tree from the word “knickknack.” The overall word with its associated indexes is shown in (A). (B) shows the unsorted suffix array with their associated indexes. (C) shows the completed suffix tree.

Figure 3

Alignment rate distributions for each of the seven aligners tested across 48 RNA-seq samples. Aligners with the same character above the boxplots were not significantly different based on Tukey’s HSD with 95% confidence level.

Figure 4

Runtime per read distributions for each of the seven aligners tested across 48 RNA-seq samples. Aligners with the same character above the boxplots were not significantly different based on Tukey’s HSD with 95% confidence level.

Figure 5

Transcriptome coverage based on a minimum alignment length cutoff. This figure shows the transcriptome coverage for each aligner calculated from the BLAST+ alignment results. The y-axis shows the transcriptome coverage calculated for varying alignment length cutoffs ranging from 100 to 2,500 base pairs (x-axis). Coverage for an aligner was calculated taking the total number of alignments returned by BLAST+ with length greater than or equal to the cutoff value and dividing by the total number of transcripts in the Erysiphe necator reference transcriptome. Each line is a different aligner except for the “Reference” plot, which represents the percent of transcripts from the reference transcriptome that follow the alignment length cutoffs. This line represents the theoretical maximum for each cutoff value.

Figure 6

Average runtime per read vs. average alignment rate for modern aligners. TopHat2 was also used but the results were excluded from this figure as it was found to be an outlier. The x-axis shows the average percentage of reads successfully mapped to the E. necator reference genome from all 48 samples. The y-axis shows the average runtime per read from all 48 samples.