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 Seq Gene Expression: A Practical Guide

RNA sequencing (RNA seq) gene expression is the process of using high-throughput sequencing to measure the abundance of RNA transcripts in a biological sample. This guide explains the core concepts, decision points, workflow, quality checks, common mistakes, and interpretation limits of RNA seq gene expression analysis. It is intended for researchers, bioinformatics beginners, and lab scientists who need a clear, source bounded framework to design, execute, and evaluate RNA seq experiments. For a foundational overview, the EMBL EBI Training materials provide excellent background on the technology. If you are running your first analysis, the Galaxy Training Network offers hands on tutorials that complement this guide.

RNA seq has replaced microarrays as the standard method for transcriptome wide expression profiling. It captures both known and novel transcripts, detects splice variants, and provides a digital count of gene expression levels. This guide covers bulk RNA seq analysis. Single cell RNA seq follows a different workflow and is referenced separately.

At a Glance

Aspect Description
Purpose Quantify gene expression levels from RNA transcripts
Typical data FASTQ files from high throughput sequencing (Illumina most common)
Key steps Quality control, trimming, alignment, quantification, normalization, differential expression, functional enrichment
Main tools FastQC, Cutadapt, STAR or HISAT2, featureCounts or Salmon, DESeq2, edgeR
Output Gene count matrix, list of differentially expressed genes, enrichment results
Requirements Reference genome or transcriptome, at least 20 million reads per sample (varies)
Typical runtime Hours to days depending on sample size and computing resources

The Bioconductor project provides validated R packages for most steps. The NCBI Sequence Read Archive is the primary public repository for raw RNA seq data.

Decision Criteria

Use RNA seq for gene expression when you need to discover new transcripts, quantify expression across the entire transcriptome, or work with organisms lacking commercial microarrays. Consider alternative methods such as qPCR for targeted validation or microarrays for very large sample sets when budget is tight. RNA seq is also superior for detecting low abundance transcripts and alternative splicing.

Choose stranded library preparation if you need to know which strand the transcript originates from. Unstranded libraries are cheaper but lose this information. Paired end reads give higher mapping accuracy and better splice junction detection compared to single end reads. Poly A selection enriches for mRNA and reduces ribosomal RNA contamination. Ribodepletion is better if you want to capture non coding RNAs.

For differential expression analysis, the Cell Mol Biol Lett study used a design with at least three biological replicates per condition. The number of replicates directly affects statistical power. Aim for a minimum of three, but more replicates (six or more) improve detection of small fold changes.

Practical Workflow and Implementation Steps

Step 1: Raw Data and Quality Control

Obtain FASTQ files from your sequencing facility or download from NCBI Sequence Read Archive. Run FastQC on each file to assess base quality scores, GC content, adapter contamination, and sequence duplication levels. Based on the FastQC report, trim or remove low quality bases and adapters using Cutadapt or Trimmomatic. Always document trimming parameters. The Galaxy Training Network provides a complete quality control tutorial.

Step 2: Alignment to Reference

Align trimmed reads to the reference genome or transcriptome with a spliced aligner such as STAR or HISAT2. STAR is faster for larger genomes but uses more memory. HISAT2 works well with limited RAM. For organisms with a well annotated reference genome, align to the genome. For non model organisms or when you lack a genome, align to a de novo assembled transcriptome using Salmon or Kallisto in quasi mapping mode. These alignment free tools are also faster and produce count data directly.

Quality check the alignment: the percentage of uniquely mapped reads should exceed 70% for most samples. Low mapping rates indicate contamination, poor library quality, or mismatched reference. Use SAMtools flagstat to verify mapping statistics. The EMBL EBI Training covers alignment evaluation.

Step 3: Quantification

Count reads overlapping gene features. Use featureCounts from the subread package or HTSeq count. If you used Salmon or Kallisto, you already have transcript level counts. Import transcript level counts to gene level using tximport in R. The output is a count matrix with rows as genes and columns as samples.

Step 4: Normalization and Differential Expression

Normalize counts to account for sequencing depth and library composition. Common methods include TPM (transcripts per million) for within sample comparison, and DESeq2 median of ratios or edgeR trimmed mean of M values (TMM) for between sample comparisons. Use DESeq2 or edgeR in R for differential expression analysis. Design a model matrix with relevant covariates (e.g., condition, batch). Filter lowly expressed genes before running.

Interpret results using adjusted p values (Benjamini Hochberg) and log2 fold changes. A typical threshold is adjusted p value < 0.05 and absolute log2 fold change > 1 (2 fold change). The Naunyn Schmiedebergs Arch Pharmacol study used these criteria to identify relevant genes.

Step 5: Functional Enrichment and Interpretation

Use tools like clusterProfiler, DAVID, or GSEA to test whether differentially expressed genes are enriched in specific Gene Ontology terms, KEGG pathways, or other gene sets. Visualize results with volcano plots, heatmaps, and MA plots. The Int J Biol Macromol study showed how chromatin interaction analysis can complement expression data. The Discov Oncol study linked RNA binding protein targets to expression changes in colorectal cancer.

Quality Checks Summary

  • Raw data: base quality, adapter content, duplication levels.
  • Trimming: no over trimming, verify loss of reads is less than 10%.
  • Alignment: uniquely mapped reads > 70%, low multimapping.
  • Quantification: high correlation between replicate samples (Pearson R > 0.9).
  • Differential expression: p value distribution is uniform for non DE genes, MA plot shows no low count bias.
  • Enrichment: adjust for multiple testing, validate top hits with qPCR if possible.

Common Mistakes

  1. Insufficient Replicates: Using only two replicates per condition yields unreliable results. Three is the minimum but more is better.
  2. Ignoring Batch Effects: If samples are processed in different batches, include batch as a covariate in your model. Use principal component analysis to check for batch clustering.
  3. Over trimming: Aggressive trimming can remove biological content. Always check that trimmed reads still have acceptable length.
  4. Choosing the Wrong Normalization: Do not use RPKM or FPKM for differential expression. Use DESeq2 or edgeR recommended normalization.
  5. Misinterpreting Fold Changes: A small fold change with very low expression may not be biologically meaningful. Filter based on expression level.
  6. Using Default Parameters Without Checking: Software defaults may not suit your data. Adjust parameters based on library type (stranded vs unstranded), read length, and insert size.
  7. Not Documenting the Analysis: Reproducibility requires recording exact software versions, commands, and parameter settings. Use a workflow manager or script.

The NCBI Bookshelf offers detailed discussions on experimental design pitfalls.

Limits and Interpretation

RNA seq measures relative expression, not absolute transcript number. Differences in expression between genes within the same sample are not directly comparable because of transcript length biases. Normalization methods correct for this but do not make expression levels absolute.

Counts are affected by technical factors: library preparation efficiency, PCR duplication rate, and sequencing depth. Lowly expressed genes are more variable. The detection limit depends on depth. For a typical 30 million read library, genes with fewer than 10 counts are often considered unreliable.

RNA seq does not directly measure protein abundance. Post transcriptional regulation can cause discordance between RNA and protein levels. Corroborate key findings with orthogonal methods such as qPCR or western blot.

Batch effects can obscure biological signals. If samples cannot be processed simultaneously, use a balanced block design and include a batch term in your model. For retrospective data, use ComBat or limma removeBatchEffect.

When using a reference genome, discrepancies between the reference and your sample genome (e.g., from a different strain or population) can cause mapping bias. Map to the best available reference and consider using SNP tolerant aligners.

The BMC Biol study on olfactory rosettes in Senegalese sole demonstrated that tissue specific plasticity can be detected with proper experimental design. However, such findings are context dependent and may not generalize. The Mol Cancer Ther study on melanoma showed that epigenetic silencing can be reversed by decitabine, but only a subset of genes responded. Always interpret differential expression in the context of the biological system and treatment.

Frequently Asked Questions

Q1: What is the minimum sequencing depth for gene expression analysis? For standard differential expression of protein coding genes, 20 to 30 million reads per sample is typical. Lower depth (10 million) may suffice for abundant transcripts but reduces power for lowly expressed genes. Deeper sequencing (50 million) is needed for splice variant detection or rare transcripts.

Q2: Can I use RNA seq for non model organisms without a reference genome? Yes, you can perform de novo transcriptome assembly using tools like Trinity or rnaSPAdes. Then quantify expression with Salmon or Kallisto using the assembled transcriptome as reference. The quality of assembly affects accuracy.

Q3: How do I handle batch effects in RNA seq data? Include batch as a factor in your linear model (DESeq2 or edgeR). Alternatively, use the limma package with removeBatchEffect for visualization. Ensure samples from each condition are distributed across batches, not all in one batch.

Q4: What is the difference between bulk RNA seq and single cell RNA seq? Bulk RNA seq measures average expression across thousands of cells. Single cell RNA seq profiles individual cells, revealing cell type heterogeneity. Single cell data require specialized normalization (SCTransform, sctransform) and have higher dropout rates. The Cell Mol Biol Lett study used single cell and spatial transcriptomics to reveal tumor microenvironment differences.

References and Further Reading

Related Articles