MMR: A tool for Read Multi-Mapper Resolution

TL;DR

The MMR tool improves the accuracy of genome and transcriptome sequencing data analysis by resolving ambiguous read mappings.

Contribution

MMR introduces a novel method for resolving multi-mapped reads using coverage density to enhance downstream analysis accuracy.

Findings

Using MMR increases transcript quantification accuracy by 17% for reads of length 51.

MMR reduces alignment file sizes by over 50%, decreasing quantification tool runtime.

The MMR algorithm is efficient, linear in complexity, and compatible with BAM files.

Abstract

Mapping high throughput sequencing data to a reference genome is an essential step for most analysis pipelines aiming at the computational analysis of genome and transcriptome sequencing data. Breaking ties between equally well mapping locations poses a severe problem not only during the alignment phase, but also has significant impact on the results of downstream analyses. We present the multimapper resolution (MMR) tool that infers optimal mapping locations from the coverage density of other mapped reads.

Filtering alignments with MMR can significantly improve the performance of downstream analyses like transcript quantitation and differential testing. We illustrate that the accuracy (Spear-man correlation) of transcript quantification increases by 17% when using reads of length 51. In addition, MMR decreases the alignment file sizes by more than 50% and this leads to a reduced running time of the quantification tool. Our efficient implementation of the MMR algorithm is easily applicable as a post-processing step to existing alignment files in BAM format. Its complexity scales linearly with the number of alignments and requires no further inputs.

Source code and documentation are available for download at github.com/ratschlab/mmr. Supplementary text and figures, comprehensive testing results and further information can be found at bioweb.me/mmr.

akahles@cbio.mskcc.org and raetsch@cbio.mskcc.org

Figures

MMR Results on Simulated Read Data:A MMR effect on artificial A. thaliana genome tiling data. Distribution of unfiltered coverage values (left) and MMR-filtered values (right). B Quantification results on 3 million simulated 51nt reads aligned with PALMapper, using rQuant and Cufflinks. Accuracy was measured as Spearman correlation to ground truth. Unfiltered read sets are shown in dark, best-hit read sets in medium and MMR -filtered in light gray. Native error rate is 2.7%, additional noise levels 1%, 2% & 3%.

Equations

Full text

Introduction

1

Addressing the increasing need for fast and accurate mapping of high throughput sequencing data to a reference sequence, many different software tools have been developed over the past years, many of which are frequently updated and improved [10, 6, 3, 7]. While numerous challenges have been addressed by the developers, e.g., the consideration of gaps and mismatches or the spliced alignment of RNA-Sequencing data, the problem of ambiguous read mapping still remains unresolved for many of the most popular alignment tools. Depending on factors like read length, alignment sensitivity, and repetitiveness of the target genome, a large fraction of reads aligns uniquely to the target and exactly one mapping location is reported. However, for the remaining, still significantly large, fraction of reads (*≈*10−20%, depending on alignment sensitivity), several possible mapping locations exist. Currently, different strategies are employed to deal with these reads in downstream analyses, most of which have unfavorable side effects: Discarding reads with ambiguous alignments from the alignment result leads to a systematic underestimation of abundance in genomic regions with multi-mapper ambiguities, whereas picking a random alignment or distributing weight across all alignments uniformly does not have a proper biological justification.

Here, we present a simple, yet powerful tool, called the Multi-Mapper Resolution tool (MMR), that assigns each read to a unique mapping location in a way that the overall read coverage across the genome is as uniform as possible. MMR makes use of the critical fraction of unambiguously aligned reads and iteratively selects the alignments of ambiguously mapping reads in a way the overall coverage becomes more uniform. We show that this strategy has a positive influence on downstream analyses, such as transcript quantification and prediction.

Approach

2

Outline of Algorithm

2.1

Our approach to resolve ambiguous mapping locations is based on the simple assumption that, besides all existing biases from library preparation and sequencing, the alignment coverage should generally be uniform within a region (RNA-seq or whole exome-seq) or the whole genome (WGS-seq). Based on this assumption, we can evaluate the fit of an alignment of a read to its current mapping location relative to other locations, by assessing the local coverage of the candidate regions. For each read the algorithm jointly evaluates all available alignments with the goal of selecting the alignment/mapping that results in the smoothest overall coverage. At the beginning, for each read one alignment is selected based on either best alignment score, the given input order or random choice. The set of all initially picked alignments as well as alignments of uniquely mapped reads define a global coverage map. Based on this map, we can evaluate the quality of an alignment in its coverage context. To choose the locally optimal alignment for each read, we perform a comparison of all alignments a with respect to a loss function ℓ^+^(a) of placing a relative to not placing it (ℓ^−^(a)). In the simplest case the loss function is defined as the empirical variance of the read coverage within window around the alignment (see Suppl. Material). This quantity can be computed efficiently since we keep track of the global coverage map, which is updated when the selected alignment changes. Given the currently selected alignment a and an alternative alignment b, we update our choice, if the overall loss across the genome would be reduced by choosing the alternative alignment. This is the case when ℓ^−^(a) + ℓ^+^(b) < ℓ^+^(a) + ℓ^−^(b). This is repeated for all reads with ambiguous mapping locations. Several iterations over the whole alignment file improve the results. However, the most striking improvements are achieved within the first three iterations and only slight changes can be observed after that. A more detailed description is provided in the Suppl. Methods section.

Paired-end Reads

2.2

Handling paired-end reads in our framework is straightforward: Instead of evaluating two individual mapping locations, the same principle is used to compare two pairs of alignments. After generating a list of proper pairs, where feasibility is determined through orientation, chromosome and reasonable proximity, the list is evaluated the same way as the list of possible mapping locations for single-end reads. This approach is easily adaptable to n-tuple of mates for n > 2.

Limiting Ambiguity

2.3

To find a good trade-off between mapping sensitivity and the number of possible mapping locations, we allow to restrict the list of possible mapping locations. This is achieved by thresholding the difference in edit operations between the best hit and any other alignment. For instance, a filter of 0 would only include alignments as possible mapping locations that have as few edit operations as the best mapping.

Implementation

2.4

The MMR approach is implemented in C++ and its source code is publicly available for download under http://github.com/ratschlab/mmr. Although it has been tested and optimized for Linux based systems, it can be compiled on other platforms. Parsing of alignment files in BAM format requires samtools [9]. We also provide a multi-threaded implementation that keeps the coverage information in common memory, requiring no additional memory if multiple threads are used. The single threaded running time depends on the number of possible mapping locations per read but is on average 30−45 seconds per one million alignments per iteration. Thus, running MMR for three iterations on 100 million alignments takes *≈*20 minutes using 10 threads (Intel Xeon E5-2665 CPU).

Application

3

As a proof of principle, we tiled the A. thaliana genome with overlapping 50-mers and aligned these 50nt reads back to the genome. This resulted in a non-uniform coverage, in particular near repetitive regions (see Fig. 1A; Suppl. Fig. 3). Using MMR, we could fully resolve mapping ambiguities, resulting in the expected uniform coverage of 50 almost everywhere.

MMR Results on Simulated Read Data:A MMR effect on artificial A. thaliana genome tiling data. Distribution of unfiltered coverage values (left) and MMR-filtered values (right). B Quantification results on 3 million simulated 51nt reads aligned with PALMapper, using rQuant and Cufflinks. Accuracy was measured as Spearman correlation to ground truth. Unfiltered read sets are shown in dark, best-hit read sets in medium and MMR -filtered in light gray. Native error rate is 2.7%, additional noise levels 1%, 2% & 3%.

We further tested our approach on a set of 3 million artificial RNA-Seq reads that were generated with FluxSimulator [5] based on a subset of 5, 000 genes randomly drawn from the human ENSEMBL annotation. We simulated 51nt and 76nt reads, resulting in an average coverage of 18*×* and 25*×*, respectively. The reads were then mutated using an empirical error model that led to a native error rate of 2.7%. Three levels of random noise (+1%, +2, +3%) were applied in addition. We aligned the reads with TopHat2 (v2.0.2; [7]) and PALMapper (v0.5; [6]), allowing for up to 6 edit operations, with no annotation provided. Further information is provided in the Suppl. Material. To investigate the effect of MMR on downstream analyses, we performed transcript quantification using Cufflinks [12] (v1.3) and rQuant [2] on the MMR-filtered alignments, the best alignments only (the alignment ranked highest by aligner) and on completely unfiltered alignments. For TopHat2 as well as PALMapper the quantifications based on the MMR-filtered alignments showed a consistently better correlation to the ground truth than both the best-hit and unfiltered alignments sets. The shorter reads of length 51nt (Fig. 1B) showed larger improvements compared to unfiltered (Cufflinks : 14.8%, rQuant : 16.6%) and best-hit set (Cufflinks : 3.2%, rQuant : 3.0%) than the longer reads of length 76nt, that showed consistent but smaller improvements (Suppl. Fig. 4 & 5).

Conclusion

4

We presented MMR, a post-processor for BAM files, resolving ambiguous alignment locations. We showed its easy applicability to the output of different alignment methods and illustrated that MMR can greatly improve accuracy of downstream quantification methods. Whereas the improvements seem moderate on a global scale, the effect on single genes can be much larger. Given its lean implementation and the short running time, MMR is very well suited for large-scale genome-, exome- and RNA-sequencing efforts. It may also be useful for post-processing alignments in meta-genome projects for improved selection of taxa.