DNA-protein quasi-mapping for rapid differential gene expression analysis in non-model organisms

Abstract

Background: Conventional differential gene expression analysis pipelines for non-model organisms require computationally expensive transcriptome assembly. We recently proposed an alternative strategy of directly aligning RNA-seq reads to a protein database, and demonstrated drastic improvements in speed, memory usage, and accuracy in identifying differentially expressed genes.

Result: Here we report a further speed-up by replacing DNA-protein alignment by quasi-mapping, making our pipeline > 1000× faster than assembly-based approach, and still more accurate. We also compare quasi-mapping to other mapping techniques, and show that it is faster but at the cost of sensitivity.

Conclusion: We provide a quick-and-dirty differential gene expression analysis pipeline for non-model organisms without a reference transcriptome, which directly quasi-maps RNA-seq reads to a reference protein database, avoiding computationally expensive transcriptome assembly.

Keywords: DNA-protein alignment; Differential gene expression analysis; Non-model organism; Quasi-mapping; RNA-seq.

Fig. 1

Example of the reference index data-structures for the set of sequences $\mathtt{MVVAV},\mathtt{VAVNV},\mathtt{VAAVV}$. C is the concatenated string, B is the bit-vector, SA suffix array, and HT is a hash table mapping a 2-mer to a suffix array interval containing suffixes whose prefix is the 2-mer

Algorithm 1

Quasi-mapping of one translation of a read. This process is repeated for all possible 6 frames.

Fig. 2

Mapping performance of the different aligners/mappers when using different reference proteomes. For DIAMOND, different points correspond to the different available presets: normal, fast, mid-sensitive, sensitive, more-sensitive, very-sensitive, and ultra-sensitive. For Kaiju, the points correspond to setting the mode to greedy and MEM. For LAST, the maximum mismap probability was set to 0.95,0.9, and 0.8. For our method, we used the k-mer size of 7 and coverage threshold was varied among 40, 50, and 60

Fig. 3

Comparison of mapping run times for a typical sample containing roughly 21 million pairs reads of length 100 bp each. Each tool was run on a single thread and with the following settings: normal for DIAMOND, greedy for Kaiju, maximum mismap probability of 0.95 for LAST, and k-mer 7 with coverage threshold 40 for our method

Fig. 4

Precision-recall curves for differential gene expression analysis using different tools for the mapping step and when using different reference proteomes. The three shape markers in each curve correspond to setting the false discovery rates in DESeq2 to the values of 0.01, 0.05, and 0.1

Fig. 5

Comparing the performance of our method on the full amino acid alphabet versus a reduced one of size 11. For each curve, the three points from left to right correspond to coverage thresholds of 60, 50, and 40, respectively