Assembly-free rapid differential gene expression analysis in non-model organisms using DNA-protein alignment

Abstract

Background: RNA-seq is being increasingly adopted for gene expression studies in a panoply of non-model organisms, with applications spanning the fields of agriculture, aquaculture, ecology, and environment. For organisms that lack a well-annotated reference genome or transcriptome, a conventional RNA-seq data analysis workflow requires constructing a de-novo transcriptome assembly and annotating it against a high-confidence protein database. The assembly serves as a reference for read mapping, and the annotation is necessary for functional analysis of genes found to be differentially expressed. However, assembly is computationally expensive. It is also prone to errors that impact expression analysis, especially since sequencing depth is typically much lower for expression studies than for transcript discovery.

Results: We propose a shortcut, in which we obtain counts for differential expression analysis by directly aligning RNA-seq reads to the high-confidence proteome that would have been otherwise used for annotation. By avoiding assembly, we drastically cut down computational costs - the running time on a typical dataset improves from the order of tens of hours to under half an hour, and the memory requirement is reduced from the order of tens of Gbytes to tens of Mbytes. We show through experiments on simulated and real data that our pipeline not only reduces computational costs, but has higher sensitivity and precision than a typical assembly-based pipeline. A Snakemake implementation of our workflow is available at: https://bitbucket.org/project_samar/samar .

Conclusions: The flip side of RNA-seq becoming accessible to even modestly resourced labs has been that the time, labor, and infrastructure cost of bioinformatics analysis has become a bottleneck. Assembly is one such resource-hungry process, and we show here that it can be avoided for quick and easy, yet more sensitive and precise, differential gene expression analysis in non-model organisms.

Keywords: DNA-protein alignment; Differential gene expression analysis; Non-model organisms; RNA-seq; Transcriptome assembly.

Fig. 1

(Left) Precision-recall curves of our method using the D. melanogaster proteome reference and the Bowtie2-RSEM-DESeq2 pipeline using its transcriptome reference. The three shape markers in each curve correspond to setting the FDR threshold of DESeq2 to 0.01, 0.05, and 0.1. (Right) Log fold change of true positive DE genes estimated by DESeq2 at FDR threshold of 0.1, compared against the 6 simulated log-fold change levels

Fig. 2

Precision-Recall curves of our method and the assembly-based pipeline, when using reference proteomes of close relatives (D. ananassae and D. grimshawi) and a distant relative (A. gambiae). The three shape markers in each curve correspond to setting the FDR values of 0.01, 0.05, and 0.1

Fig. 3

Precision-Recall curves of our method and the assembly-based pipeline, when using reference proteomes of close relatives (D. ananassae and D. grimshawi) and a distant relative (A. gambiae), and with the alignments produced by Dammit which covered less than 50% of the length of the contig removed. The three shape markers in each curve correspond to setting the FDR values of 0.01, 0.05, and 0.1

Fig. 4

For the three reference proteomes, Venn diagrams showing the intersections among the Baseline set consisting of DE genes called by the baseline approach of Bowtie2-RSEM-DESeq2 using the D. melanogaster reference transcriptome, Our set consisting of D. melanogaster genes to which the DE genes called by our approach mapped to, and (3) Assembly-based set consisting of D. melanogaster genes to which Trinity DE genes mapped to. FDR threshold of 0.01 was used for all three approaches