Zubair Khalid

Virologist/Molecular Biologist | Veterinarian | Bioinformatician

Conventional & Molecular Virology • Vaccine Development • Computational Biology

Dr. Zubair Khalid is a veterinarian and virologist specializing in conventional and molecular virology, vaccine development, and computational biology. Dedicated to advancing animal health through innovative research and multi-omics approaches.

Dr. Zubair Khalid - Veterinarian, Virologist, and Vaccine Development Researcher specializing in Computational Biology, Multi-omics, Animal Health, and Infectious Disease Research

Category: Guides

Rna Sequencing Analysis Pipeline

An RNA sequencing analysis pipeline is a standardized sequence of computational steps that transform raw next generation sequencing reads into quantifiable transcriptome data, typically gene expression counts or transcript abundance. This guide is intended for research scientists, bioinformatics students, and laboratory staff who need a practical, source bounded framework to design, execute, and interpret an RNA seq analysis without relying on opaque or all purpose workflows. Every major decision point, quality check, and interpretative limit is grounded in authoritative technical references such as the NCBI Bookshelf and EMBL EBI Training resources, ensuring that your results are reproducible and biologically meaningful.

For example, the NCBI Bookshelf provides foundational reviews of sequencing technologies and quality metrics [NCBI Bookshelf]. Meanwhile, the EMBL EBI Training portal offers detailed tutorials on read alignment and differential expression analysis that align with community best practices [EMBL EBI Training]. This guide builds on those resources to deliver a concrete, step by step approach.

At a Glance

Step Purpose Representative Tool Essential Quality Check
Raw data acquisition Obtain reads from public repositories or local runs NCBI SRA, local FASTQ Verify read count, base quality scores, adapter contamination
Quality control (QC) Assess read quality and detect problems FastQC Per base sequence quality, GC content, overrepresented sequences
Read trimming Remove adapters and low quality bases Trimmomatic, Cutadapt Compare pre and post trimming QC reports
Alignment or pseudoalignment Map reads to a reference genome or transcriptome STAR, HISAT2, Salmon Alignment rate, uniquely mapped percentage, uniformity of coverage
Quantification Count reads per gene or transcript featureCounts, Salmon quant Number of detected features, zero count distribution
Normalization Adjust for library size and composition DESeq2, edgeR, limma voom Library size distribution, PCA or MDS plot
Differential expression Identify significantly changed genes DESeq2, edgeR Dispersion estimates, p value distribution, fold change distribution
Functional interpretation Enrichment analysis, pathway mapping clusterProfiler, GSEA Gene ratio, adjusted p value, biological theme overlap

Decision Criteria

Every RNA seq experiment requires explicit choices at the pipeline design stage. The most critical decisions include the strandedness of the library, the read length and paired end status, the choice between genome based alignment and transcriptome based pseudoalignment, and the specific quantification approach.

Strandedness determines which strand of the original transcript is captured. Stranded protocols preserve the orientation of reads relative to the transcript, while unstranded libraries lose that information. For differential expression analysis of protein coding genes, unstranded data is often sufficient, but for antisense transcripts, noncoding RNAs, or isoform level quantification, stranded data substantially improves accuracy. The Galaxy Training Network provides clear examples of how strandedness affects alignment counts [Galaxy Training Network].

Read length and paired end sequencing affect splice junction detection. Longer reads (at least 75 bp) improve junction mapping, and paired end reads increase mapping confidence across splice sites. For standard differential expression of annotated genes, single end 50 bp reads can still be adequate, but for novel isoform discovery or low abundance transcripts, paired end 100 bp or longer is recommended. The Bioconductor project includes vignettes that benchmark these parameters for real RNA seq experiments [Bioconductor].

The choice between a genome based aligner like STAR and a pseudoaligner like Salmon often hinges on the experimental goal. Genome based alignment is necessary for detecting novel splice junctions or fusion transcripts. Pseudoalignment is faster and requires less memory, making it popular for quantifying known transcripts in large studies. However, pseudoalignment assumes a complete and accurate transcriptome annotation. If your experiment involves a nonmodel organism or an incomplete annotation, genome based alignment with featureCounts is more reliable. The EMBL EBI Training materials offer side by side comparisons of these strategies.

Practical Workflow and Implementation Steps

Implement a robust RNA seq pipeline by following these ordered steps. Each step includes a built in quality control to catch errors early.

  1. Obtain raw data. Download FASTQ files from the NCBI Sequence Read Archive or generate them from your own sequencing run. Verify the number of reads and base quality immediately [NCBI SRA].

  2. Quality control. Run FastQC on all raw samples. Look for per base sequence quality, adapter contamination, and GC bias. Flag any sample with a high proportion of reads below Q20.

  3. Read trimming. Use Trimmomatic to remove adapter sequences and trim low quality bases from the ends. Use a sliding window approach (e.g., window size 4, average quality 15). After trimming, rerun FastQC to confirm improvement.

  4. Alignment to reference. Choose a reference genome or transcriptome from Ensembl or RefSeq. For genome alignment, use STAR with options for outFilterMultimapNmax, outFilterMismatchNmax, and twopassMode if detecting novel transcripts. Assess alignment rate (typically > 80% for good data). If using Salmon for pseudoalignment, provide an index built from the transcriptome FASTA. Check that the mapping rate matches library complexity.

  5. Gene level quantification. For STAR aligned reads, run featureCounts to assign reads to genes using a GTF annotation file. Ensure that the strand specific option is set correctly (0 for unstranded, 1 for stranded, 2 for reversely stranded). For Salmon, the quant.sf files contain transcript level counts that can be aggregated to gene level using tximport in R.

  6. Normalization and exploratory analysis. Import count data into R using DESeq2. Use DESeq2's median of ratios normalization to account for library size and composition differences. Generate a PCA plot to check for batch effects or sample outliers. If a batch effect is visible, include batch as a covariate in the design formula.

  7. Differential expression analysis. Run DESeq2 to compute log2 fold changes and adjusted p values (Benjamini Hochberg). Filter for genes with an adjusted p value <= 0.05 and a fold change above a biologically meaningful threshold (e.g., |log2FC| > 1). Check the dispersion plot and p value histogram for anomalies.

  8. Functional interpretation. Use clusterProfiler or GSEA to perform gene ontology enrichment or KEGG pathway analysis. Consider the background of expressed genes to avoid bias. Report only pathways that pass a multiple testing correction.

Each of these steps can be validated by comparing with the official Galaxy Training Network RNA seq tutorials, which provide hands on exercises for human and mouse data [Galaxy Training Network]. Additionally, the Bioconductor workflow for RNA seq analysis offers reproducible code examples that follow the same logic [Bioconductor].

Common Mistakes

Several recurrent errors undermine the reliability of RNA seq results. The most widespread mistake is ignoring batch effects. Even with careful library preparation, technical variation across sequencing runs can mask biological signals. Always include batch in the design matrix, and run a PCA before the final model.

Another common mistake is using the wrong strandedness flag. Forcing an unstranded alignment on stranded data produces inflated counts for intronic and antisense reads, biasing differential expression. Verify the strandedness by comparing alignment statistics to the library preparation protocol.

Researchers also frequently use a single reference annotation without checking its version. Gene symbols and transcript identifiers change between releases. Using a different annotation for alignment and quantification leads to mismatched counts. Always download the annotation from the same source and version as the reference genome.

A third mistake is ignoring multimapping reads. Reads that map to multiple locations, often from paralogous genes or repetitive regions, should be treated consistently. Many pipelines assign multimappers proportionally or discard them. The choice significantly affects expression estimates for gene families. The NCBI Bookshelf discusses multimapping strategies in detail [NCBI Bookshelf].

Finally, over interpreting small fold changes without biological validation is a common pitfall. Even with a low p value, a 1.2 fold change may not be reproducible or functionally relevant. Always perform RT qPCR or another orthogonal method on a subset of candidate genes.

Limits of Interpretation

RNA seq data provide relative, not absolute, transcript abundance. A gene with higher normalized counts in one condition is expressed more than in another condition, but the pipeline cannot determine the absolute number of transcripts per cell. Moreover, RNA seq measures steady state RNA levels, which combine transcription, processing, and degradation. Changes in expression may reflect altered RNA stability rather than transcriptional regulation.

Interpretation is also constrained by sequencing depth. Shallow sequencing (under 10 million reads per sample) may fail to detect low abundance transcripts and can inflate noise in differential expression. The Single cell RNA sequencing study of aortic dissection shows how depth impacts detection of T cell transcripts [pubmed 42436415]. Similarly, the cross species analysis of Mycobacterium cavoris required deep sequencing to identify rare transcripts [pubmed 42435264].

Another important limit: alternative isoform quantification remains uncertain, especially for isoforms that share many exons. Most pipelines assign reads to the canonical isoform by default. The paper on noncoding variants in hereditary cancer syndromes notes that overlooking minor isoforms can miss disease relevant variants [pubmed 42433780]. Similarly, the simultaneous transcriptomic analysis of host and symbiont requires careful filtering to avoid cross mapping artifacts [pubmed 42427455].

Finally, the interpretation of fold changes depends on the normalization method. The sequencing level spatial simulation tool STGBench warns that different normalization strategies can produce divergent rankings of differentially expressed genes [pubmed 42418828]. Always compare outcomes from two normalization methods (e.g., FPKM vs. DESeq2) and confirm that top hits are robust.

The limits do not invalidate RNA seq. Rather, they demand cautious biological validation and transparent reporting of pipeline parameters. The multiomic profiling of murine pancreas demonstrates that integration with epigenetic data can overcome some of these limits by providing independent evidence of regulation [pubmed 42414525].

Frequently Asked Questions

What is the minimum number of reads per sample for standard differential expression?
For most human or mouse experiments, a minimum of 10 million uniquely mapped reads per sample is recommended. Fewer reads may still detect highly expressed genes but will miss lowly expressed ones and increase the false discovery rate. The EMBL EBI Training guide includes depth versus power curves.

Should I use Unique Molecular Identifiers (UMIs) in my RNA seq pipeline?
UMIs are essential for experiments that require absolute molecule counting, such as single cell RNA seq or ultra low input bulk RNA seq. For standard bulk RNA seq with high input RNA, UMIs are not necessary because PCR duplicates are rare and can be accounted for by deduplication tools. Refer to the Bioconductor documentation for UMI aware workflows.

How do I handle reads that map to multiple locations (multimappers)?
Multimapping reads are common for gene families and repetitive elements. The safest approach is to use a pseudoaligner that resolves multimappers probabilistically (like Salmon), or to assign multimappers fractionally to each locus. Uniquely mapping only will underestimate expression of multigene families. The Galaxy Training Network module on RNA seq discusses multimapping rates.

What should I do if my alignment rate is below 70%?
Check for adapter contamination by running FastQC after trimming. If the problem persists, verify that the reference genome matches the species and strain of your sample. For samples with high rRNA carryover, remove rRNA reads computationally using Bowtie2 before alignment. The NCBI Bookshelf troubleshooting section lists these steps.

References and Further Reading

Related Articles