RNA-sequencing (RNA-seq) has a wide range of applications, and there is no optimal pipeline for all cases. We review all of the major steps in RNA-seq data analysis, including quality control, read alignment, quantification of gene and transcript levels, differential gene expression, functional profiling, and advanced analysis. They will be discussed later.
Figure 1.
The general workflow of RNA-seq analysis.
Quality control of raw reads
Quality control of RNA-seq raw reads consists of analysis of sequence quality, GC content, adaptor content, overrepresented k-mers, and duplicated reads, dedicated to detecting sequencing errors, contaminations, and PCR artifacts. Read quality decreases towards the 3’ end of reads, bases with low quality, therefore, they should be removed to improve mappability. In addition to the quality of raw data, quality control of raw reads also includes the analysis of read alignment (read uniformity and GC content), quantification (3’ bias, biotypes, and low-counts), and reproducibility (correlation, principal component analysis, and batch effects).
Table 1. The tools for quality control of RNA-seq raw reads.
Tools | Applications |
NGSQC | Quality control of raw reads generated by Illumina platforms. |
FastQC | Quality control of raw reads generated by any platforms. |
FASTX-Toolkit | Discard of low-quality reads, trim adaptor sequences, and elimination of poor quality bases. |
Trimmonmatic | |
Picard | Quality control in read alignment, including the determination of read uniformity and GC content. |
RSeQC | |
Qualimap | |
NOISeq | Provide useful plots for quality control of count data. |
EDASeq |
Read alignment
There are generally three strategies for read alignment, genome mapping, transcriptome mapping, and de novo assembly. Regardless of whether a genome or transcriptome reference is available, reads may map uniquely or be assigned to multiple positions in the reference, which are referred to as multi-mapped reads or multireads. Genomic multireads are generally due to repetitive sequences or shared domains of paralogous genes. Transcriptome multi-mapping arises more often due to gene isoforms. Therefore, transcript identification and quantification are important challenges for alternatively expressed genes. When a reference is not available, RNA-seq reads are assembled de novo using packages such as SOAPdenovo-Trans, Oases, Trans-ABySS, or Trinity. PE strand-specific and long-length reads are preferred since they are more informative. Emerging long-read technologies, such as PacBio SMRT sequencing and Nanopore sequencing, can generate full-length transcripts for most genes.
Figure
2. Three basic strategies for RNA-seq read mapping (Conesa et al. 2016). Abbreviations: GFF, General
Feature Format; GTF, gene transfer format; RSEM, RNA-seq by Expectation Maximization.
Table 2. The comparison of genome-based and de novo assembly strategies for RNA-seq analysis.
Genome-based | De novo assembly | |
Method | Alignment to a reference genome | Not using a reference genome |
Advantages |
|
|
Disadvantages | Requires high-quality reference genome |
|
Recommended depth | Approximately 10x | Beyond 30x |
Table 3. The public sources of RNA-seq data.
Transcriptomic Database |
Data Type | Website |
Gene Expression Omnibus (GEO) | Both microarray and sequencing data | https://www.ncbi.nlm.nih.gov/geo/ |
ArrayExpress | Both microarray and sequencing data | https://www.ebi.ac.uk/arrayexpress/ |
ENCODE: Encyclopedia of DNA Elements | Public ENCODE Consortium data | https://www.encodeproject.org/ |
Sequence Read Archive (SRA) | Sequencing data | https://www.ncbi.nlm.nih.gov/sra |
European Nucleotide Archive (ENA) | Sequencing data | https://www.ebi.ac.uk/ena |
DDBJ Sequence Read Archive (DRA) | Sequencing data | https://www.ddbj.nig.ac.jp/dra |
Transcript quantification
Transcript quantification can be used to estimate gene and transcript expression levels.
Table 4. The common tools for transcript quantification.
Tools | Principles and Applications |
TopHat | Using an expectation-maximization approach that estimates transcript abundances. |
Cufflinks | Designed to take advantage of PE reads, and may use GTF information to identify expressed transcripts, or can infer transcripts de novo from the mapping data alone. |
RSEM | Quantify expression from transcriptome mapping. Allocate multi-mapping reads among transcript and output within-sample normalized values corrected for sequencing biases. |
Sailfish | |
kallisto | |
NURD | Provides an efficient way of estimating transcript expression from SE reads with a low memory and computing cost. |
Figure
3. The tools for isoform expression quantification.
Differential expression testing
Differential expression testing is used to evaluate if one gene is differentially expressed in one condition compared to the other(s). Normalizing methods need to be adopted before comparing different samples. RPKM and TPM normalize away the most important factor, sequencing depth. TMM, DESeq, and UpperQuartile can ignore highly variable and/or highly expressed features. Other factors that interfere with intra-sample comparisons involve transcript length, positional biases in coverage, average fragment size, and GC content, which can be normalized by tools, such as DESeq, edgeR, baySeq, and NOISeq. Batch effects may still be present after normalization, which can be minimized by appropriate experimental design, or removed by methods such as COMBAT or ARSyN.
Table 5. The normalization tools for differential expression testing.
Package |
Read count distribution assumptions | Input | Replicates | Normalization |
DESeq | Negative binomial distribution | Raw counts | No | Library size |
edgeR | Bayesian methods for negative binomial distribution | Raw counts | Yes | Library size TMM RLE Upperquartile |
baySeq | Bayesian methods for negative binomial distribution | Raw counts | Yes | Library size Quantile TMM |
NOISeq | Non-parametric | Raw or normalized counts | No | Library size RPKM TMM Upperquartile |
Alternative splicing analysis
Alternative splicing (AS) is a posttranscriptional process which generates different transcripts from the same gene and is vital in response to environmental stimuli by producing diverse protein products. Multiple bioinformatics tools have been developed to detect AS from experimental data. The comparison of these detection tools using RNA-seq data was conducted by Ding in 2017, and the results are shown in Table 7. They have demonstrated that TopHat and its downstream tool, FineSplice, are the fastest tools, whereas PASTA is the slowest program. Furthermore, AltEventFinder can detect the highest number of junctions, and RSR detects the lowest number of junctions. Other tools, such as TopHat, are likely to detect false positive ones. Of the two tools that detect differentially spliced isoforms, rMATS is faster than rSeqDiff but detects less differentially spliced isoforms than rSeqDiff.
Table 7. Detected AS types or differentially spliced isoforms of these tools (Ding et al. 2017).
Tool | Data Source | Running Time (Minutes) | Maximum Memory (GB) |
Maximum CPU (%) | Number of SJs | Number of Differentially Spliced Isoforms |
Alt Event Finder | ENCODE | 12 | 1.364 | 100 | 30569 | N/A |
SpliceMap | ENCODE | 42 | 3.1 | 99.9 | 11882 | N/A |
FineSplice | ENCODE | 2 | 1.364 | 100 | 8577 | N/A |
RSW | N/A | N/A | N/A | N/A | N/A | N/A |
RSR | ENCODE | 24 | 3.968 | 100 | 3143 | N/A |
PASTA | ENCODE | 350 | 2.17 | 101 | 14675 | N/A |
rMATS | mouse used in RSW study | 44 | 26.536 | 274 | N/A | 17 |
SOAPsplice | ENCODE | 123 | 5.332 | 99.7 | 10381 | N/A |
SplicePie | N/A | N/A | N/A | N/A | N/A | N/A |
SplicingCompass | N/A | N/A | N/A | N/A | N/A | N/A |
TopHat | ENCODE | 1.75 | 1.364 | 100 | 9619 | N/A |
TrueSight | ENCODE | 229 | 2.914 | 571 | 12360 | N/A |
NSMAP | N/A | N/A | N/A | N/A | N/A | N/A |
rSeqDiff | mouse used in RSW study | 115 | 0.186 | 119 | N/A | 203 |
rSeqNP | N/A | N/A | N/A | N/A | N/A | N/A |
Visualization
There are many bioinformatics tools for the visualization of RNA-seq data, including genome browsers, such as ReadXplorer, UCSC browser, Integrative Genomics Viewer (IGV), Genome Maps, Savant, tools specifically designed for RNA-seq data, such as RNAseqViewer, as well as some packages for differential gene expression analysis that enable the visualization, such as DESeq2 and DEXseq in Bioconductor. Packages, such as CummeRbund and Sashimi plots, have also been developed for visualization-exclusive purposes.
Functional Profiling
The latest step in a standard transcriptomics study is generally the characterization of the molecular functions or pathways in which differentially expressed genes are involved. Gene Ontology, Bioconductor, DAVID, or Babelomics contain annotation data for most model species, which can be used for functional annotation. As for novel transcripts, protein-coding transcripts can be functionally annotated using orthology with the help of databases such as SwissProt, Pfam, and InterPro. Gene Ontology (GO) allows for some exchangeability of functional information across orthologs. Blast2GO is a popular tool that allows massive annotation of complete transcriptome against a variety of databases and controlled vocabularies. The Rfam database contains most well-characterized RNA families that can be used for functional annotation of long non-coding RNAs.
Advanced analysis
The advanced analysis of RNA-seq usually includes other RNA-seq and integration with other technologies, which is outlined in Figure 4. More information on applications of RNA-seq, please view this article Applications of RNA-Seq.
Figure 3. The advanced analysis of RNA-seq data.
Our experienced bioinformatics scientists are skilled in utilizing the advanced bioinformatics tools to deal with the numerous sequences generated by the next and third generation sequencing. We provide both sequencing and bioinformatics services for genomics, transcriptomics, epigenomics, microbial genomics, single-cell sequencing, and PacBio SMRT sequencing.
References: