Rna Sequencing Analysis Workflow
An RNA sequencing (RNA-seq) analysis workflow transforms raw sequencing reads into biological insights through a series of quality control, alignment, quantification, and differential expression steps. This guide provides a rigorous, source bounded framework for researchers, bioinformaticians, and students who plan to perform or evaluate RNA-seq experiments. Refer to official training from EMBL-EBI Training for foundational knowledge of the computational approaches described below.
The NCBI Sequence Read Archive (SRA) [5] is the primary public repository for raw sequencing data, and the Galaxy Training Network [3] offers interactive tutorials that mirror the workflow presented here. Throughout this guide, we emphasize practical decision points, quality checks, common pitfalls, and the limits of interpretation based on established methodologies and recent literature.
At a Glance
| Step | Key Tools | Typical Output |
|---|---|---|
| Raw data QC | FastQC, MultiQC | Quality report (per base quality, GC content, overrepresented sequences) |
| Read trimming | cutadapt, Trimmomatic | Clean FASTQ files with adapter removal and quality filtering |
| Alignment | STAR, HISAT2, Salmon (pseudoalignment) | SAM/BAM files or transcript abundance estimates |
| Quantification | featureCounts, htseq count, Salmon output | Count matrix (genes or transcripts) |
| Differential expression | DESeq2, edgeR, limma voom | List of differentially expressed genes with statistics |
| Downstream analysis | clusterProfiler, GSEA | Enriched pathways, GO terms, functional categories |
Decision Criteria
Before starting the workflow, you must resolve several experimental and computational choices:
Stranded vs. unstranded library preparation: Stranded protocols preserve the original RNA strand orientation. This choice affects quantification accuracy and is essential for correct interpretation of antisense transcription. Check your kit documentation. The EMBL-EBI Training materials explain how to detect strandness from alignment statistics.
Single end vs. paired end reads: Paired end reads improve alignment in repetitive regions and enable detection of structural variants and alternative splicing. Single end is cheaper and sufficient for gene level quantification in well annotated genomes.
Reference genome vs. transcriptome alignment: Genome alignment (with STAR or HISAT2) is necessary for detecting novel transcripts, splice junctions, or RNA editing events. Transcriptome alignment (with Salmon) is faster and often sufficient for differential expression of known genes.
Replicates: Biological replicates are essential for statistical power. A minimum of three replicates per condition is standard. Technical replicates are rarely needed for well controlled experiments.
Sequencing depth: For gene level differential expression, 10 to 30 million reads per sample is typical. Rare transcript detection or isoform analysis may require 50 million or more.
Practical Workflow or Implementation Sequence
1. Raw Data Quality Control
Inspect each FASTQ file with FastQC. Look for per base quality scores (Q > 30 expected), GC content distribution (should match expected genome composition), adapter contamination, and overrepresented sequences. Aggregate reports across samples using MultiQC. The Galaxy Training Network provides a step by step tutorial on interpreting FastQC results.
2. Read Trimming and Filtering
Remove adapter sequences and low quality bases. Use cutadapt to trim 3' adapters and quality filter (e.g., Q > 20). For paired end data, ensure both reads are trimmed together. Re run FastQC on trimmed reads to confirm improvement.
3. Alignment or Pseudoalignment
Choose between:
- Genome alignment with STAR: Faster than many alternatives, two pass mode recommended for splice junction discovery. Output a sorted BAM file.
- Pseudoalignment with Salmon: Uses a transcriptome index and provides quantification directly. Uses GC bias correction and sequence specific bias models.
For both methods, verify alignment rate (typically > 70% for good quality data) and check mapping statistics using samtools flagstat. The NCBI Bookshelf [1] contains detailed chapters on alignment principles.
4. Quantification
For genome alignment, generate gene level counts using featureCounts or htseq count. Specify strandness correctly (e.g., reversely stranded for many Illumina stranded kits). For Salmon, the quantification files can be imported into R using tximport.
A count matrix with rows as genes and columns as samples is the final input for differential expression.
5. Differential Expression Analysis
Use DESeq2 (recommended for small sample sizes and negative binomial distribution), edgeR, or limma voom. The Bioconductor [4] project provides comprehensive vignettes for each package.
Key steps:
- Normalize for library size (e.g., DESeq2 median of ratios).
- Estimate dispersion.
- Fit a model (e.g., condition + batch).
- Test for differentially expressed genes using Wald test or likelihood ratio test.
- Apply multiple testing correction (Benjamini Hochberg FDR).
Generate exploratory plots: PCA (to assess batch effects and clustering), MA plot, heatmap of top genes.
6. Downstream Functional Analysis
Use clusterProfiler or GSEA to test for enrichment of Gene Ontology terms, KEGG pathways, or other gene sets. The EMBL-EBI Training offers modules on pathway analysis interpretation.
Quality Checks at Each Stage
- After trimming: Confirm removal of adapter sequences and improvement in quality scores.
- After alignment: Check uniquely mapped percentage (should exceed 60%), mitochondrial read fraction (high may indicate cell stress or sample degradation), and junction read counts.
- After quantification: Examine library size distribution, verify that no sample has excessively low counts.
- Before differential expression: Plot sample to sample correlation. Remove outliers if justified by technical notes. Check PCA for batch effects.
- After differential expression: Inspect p value histogram (should be uniform for null, with a peak near zero). Verify that top genes show consistent expression patterns.
Common Mistakes
Ignoring strandness: Using the wrong strandness parameter in featureCounts or htseq count will produce erroneous counts. Check with a small subset or use a stranded detection script from Bioconductor.
Using default trimming without inspection: Over trimming removes bases that overlap with transcript ends. Validate with FastQC.
Incorrect replicate handling: Treating biological replicates as technical replicates inflates significance. Use proper design formula in DESeq2.
Not correcting for batch effects: If samples were processed in separate batches, include batch as a covariate in the model. Confounding batch with condition prevents correction.
Overinterpreting lowly expressed genes: Genes with very low counts have high technical noise. Filter with a minimum count threshold (e.g., 10 total counts across samples).
Forgetting to use multiple testing correction: Reporting raw p values leads to excessive false positives.
Limits of Interpretation
RNA seq differential expression reveals relative changes, not absolute transcript copy numbers. Normalization methods (e.g., TPM, FPKM) are useful for comparison within samples but are not equivalent to molecule counts.
Gene level analysis cannot detect alternative splicing. For splicing studies, use transcript level quantification with Salmon or tools like rMATS.
Results depend heavily on the quality of the reference annotation. Novel transcripts or poorly annotated genomes require de novo assembly or alignment to the genome followed by transcript discovery.
Batch effects that are confounded with the experimental condition cannot be removed statistically. Careful experimental design is essential.
The workflow described here assumes standard bulk RNA seq. Specialized applications such as single cell RNA seq [10], circulating tumor cell analysis [11], or viral amplicon sequencing [9] require modifications. For instance, long read based RNA seq [7] has distinct alignment and quantification challenges.
Recent developments include methods for detecting RNA editing sites [6] that require careful calibration to avoid false positives from sequencing errors or alignment artifacts. Similarly, AAV vector transcript characterization [8] demands specific quality controls for vector derived transcripts.
Frequently Asked Questions
What is the best aligner for RNA seq?
There is no single best aligner. STAR is widely used for speed and sensitivity with splice junctions. HISAT2 is memory efficient. Salmon uses pseudoalignment for faster quantification of known transcripts. Choose based on your reference availability and goals. The Galaxy Training Network provides benchmarks for different tools.
How do I know if my library is stranded?
Check the alignment of reads to genes and the orientation of read pairs. A simple test: align a subset to a housekeeping gene. In a stranded library, the majority of reads will align to the expected strand. Tools like infer experiment in RSeQC provide automatic detection. The EMBL-EBI Training includes a module on library quality control.
How many biological replicates are necessary?
Statistical power increases with replicates. Three replicates per condition is the minimum for reliable dispersion estimates in DESeq2. For detecting small fold changes (e.g., 1.5 fold), more replicates are needed. The Bioconductor [4] vignettes include power analysis guidance.
What is the minimum sequencing depth for differential expression?
For well expressed genes, 10 million reads per sample is often sufficient. For lowly expressed transcripts or rare isoforms, 30 to 50 million reads may be needed. Deep sequencing does not compensate for poor biological replication. Check your specific experimental context using the NCBI SRA for similar studies.
References and Further Reading
- NCBI Bookshelf , Free textbooks on sequencing technology and bioinformatics basics.
- EMBL-EBI Training , Official courses on RNA seq analysis and functional interpretation.
- Galaxy Training Network , Interactive tutorials covering each step of the workflow.
- Bioconductor , R packages for differential expression and statistical analysis.
- NCBI Sequence Read Archive , Public repository for raw sequencing data.
- High precision Detection of RNA Editing Sites using Calibrated Differential RNA Editing Scanner. PubMed 42441715 , Specialized method for RNA editing detection.
- A systematic benchmark of bioinformatics methods for single cell and spatial RNA seq nanopore long reads data. PubMed 42441075 , Evaluation of long read analysis pipelines.
- AAV vector production in suspension cells using PEI transfection and sodium butyrate with orthogonal assessment of function and quality. PubMed 42438474 , Quality control approaches for vector derived transcripts.
- A highly sensitive amplicon sequencing workflow for genomic surveillance of Usutu virus. PubMed 42426896 , Amplicon based RNA seq protocol.
- Single Cell RNA Sequencing of Pancreatic Circulating Tumor Cells. PubMed 42426462 , Application to rare cell populations.
- Characterization of Circulating Tumor Cells and Other Rare Cells in Blood Using the RareCyte Platform. PubMed 42426455 , Methodology for low input RNA seq.