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

Section: Infrastructure, Cloud & Policy

RNA Sequencing Data Analysis: From Raw Reads to Differential Expression

RNA sequencing has become a standard method for measuring gene expression across entire transcriptomes, but the journey from raw sequencing reads to biologically meaningful lists of differentially expressed genes requires a structured analytical workflow. This article outlines the core steps in RNA-seq data analysis, from quality control through alignment, quantification, and differential expression testing, with practical guidance for troubleshooting common issues such as batch effects and low mapping rates. The intended readers are students, researchers, analysts, and life-science professionals who need a reproducible framework for processing RNA-seq data and interpreting the results within the limits of the technology.

The Scope of RNA-Seq Analysis

RNA-seq analysis transforms millions of short sequencing reads into interpretable measurements of transcript abundance. The workflow begins with raw FASTQ files produced by sequencing instruments and ends with tables of genes or transcripts that show statistically significant expression changes between experimental conditions. Between these endpoints lie several critical decision points where analyst choices materially affect the final results.

A standard RNA-seq pipeline includes quality control of raw reads, preprocessing steps such as adapter trimming, alignment to a reference genome or transcriptome, quantification of gene or transcript abundance, normalization to make samples comparable, and statistical testing for differential expression. Each step has multiple software options and parameter choices, and no single pipeline works optimally for every experimental design or organism [6]. The practical implication is that analysts must understand the rationale behind each step instead of blindly following a fixed protocol.

The field continues to evolve rapidly. New analysis methods and tools appear regularly, and established workflows are frequently updated to incorporate improved algorithms [7]. Analysts should therefore treat published pipelines as starting points that require adaptation to their specific data types and biological questions.

At a Glance: Core Workflow Components

The table below summarizes the primary stages of a typical RNA-seq analysis pipeline, the common tools used at each stage, the key outputs produced, and the main quality concerns that require attention.

Workflow Stage Common Tools Primary Output Key Quality Concerns
Raw data quality control FastQC, MultiQC Per-sample quality reports Adapter contamination, low base quality scores, GC bias, overrepresented sequences
Read preprocessing Trimmomatic, cutadapt Filtered and trimmed FASTQ files Overly aggressive trimming, loss of read pairs, residual adapter sequence
Alignment STAR, HISAT2 BAM files with mapped reads Low mapping rate, high multimapping rate, alignment bias across samples
Quantification featureCounts, HTSeq-count Gene-level or transcript-level count matrices Read assignment ambiguity, multi-overlapping features, strandness errors
Normalization and differential expression DESeq2, edgeR, limma-voom Lists of differentially expressed genes Batch effects, outlier samples, library size variation, low-count filtering

The workflow stages listed here represent the conventional path for bulk RNA-seq analysis. Single-cell RNA-seq follows a related but distinct pipeline that includes cell-level quality control, dimensionality reduction, and clustering before differential expression testing [5][10]. The remainder of this article focuses primarily on bulk RNA-seq, with notes on where single-cell analysis diverges.

Experimental Design Considerations Before Analysis

Analysis decisions are constrained by experimental design choices made before sequencing begins. The number of biological replicates, the sequencing depth, the read length, and the choice of single-end versus paired-end sequencing all influence what downstream analyses are possible and how reliable the results will be.

Biological replication is the most important design factor for differential expression analysis. Technical replicates, which involve sequencing the same library multiple times, do not capture biological variation between individuals or samples. Statistical methods for differential expression require an estimate of biological variability to calculate meaningful p-values and false discovery rates. Studies with too few biological replicates will have low statistical power and may miss genuine expression differences.

Sequencing depth determines how many reads are available per sample and therefore affects the sensitivity for detecting lowly expressed genes. Deeper sequencing improves detection of rare transcripts but yields diminishing returns for highly expressed genes. The optimal depth depends on the organism, the RNA complexity, and the biological question. For differential expression of protein-coding genes in mammalian samples, tens of millions of reads per sample are commonly used, but the exact requirement varies by application.

Paired-end sequencing provides information about the physical distance between read pairs, which improves alignment accuracy for reads spanning exon junctions and enables isoform-level analysis. Single-end sequencing is less expensive and may suffice for gene-level quantification in organisms with well-annotated genomes. The choice affects downstream options for transcript discovery and splice variant analysis.

Quality Control of Raw Sequencing Reads

Quality control is the first analytical step after sequencing and serves to identify problems in the raw data before they propagate through the pipeline [7][9]. The goal is to assess whether the sequencing run produced data of sufficient quality for reliable downstream analysis and to identify any sample-level issues that require attention.

Per-Base Quality Scores

Sequencing instruments assign a quality score to each base call, representing the probability that the base is correct. FastQC and similar tools generate per-base quality plots that show the distribution of quality scores across read positions. Quality scores typically decline toward the 3-prime end of reads, which is expected for most sequencing platforms. A sudden drop in quality at a specific position may indicate an instrument problem or a batch effect.

Low-quality bases at read ends are commonly removed by trimming. The trimming threshold should balance the benefit of removing unreliable bases against the cost of shortening reads and potentially losing biological information. Analysts should examine the quality profiles across all samples in a study instead of assessing each sample in isolation, because systematic patterns may indicate run-level problems.

Adapter Contamination

Adapter sequences are added during library preparation and can appear in reads when the insert fragment is shorter than the read length. Adapter contamination is more common in libraries with short inserts and can cause alignment problems because adapter sequence does not match the reference genome. FastQC reports overrepresented sequences that often correspond to adapters. Trimming tools remove adapter sequence from reads before alignment.

The extent of adapter contamination should be recorded for each sample. High and variable adapter content across samples can introduce bias if some samples require more trimming than others. Consistent library preparation protocols minimize this variability.

Sequence Duplication and GC Bias

PCR amplification during library preparation can produce duplicate reads that represent the same original fragment. High duplication rates reduce the effective complexity of the library and can distort expression measurements. Duplicate reads are more problematic for certain applications, such as variant calling, than for standard gene-level quantification, where many analysis pipelines do not remove duplicates by default.

GC bias refers to the tendency of sequencing and amplification to favor fragments with particular GC content. This bias can affect the accuracy of expression measurements for genes with extreme GC content. The severity of GC bias varies between library preparation methods and sequencing platforms [18]. Analysts should be aware of this limitation when interpreting results for genes at the extremes of the GC distribution.

Read Alignment to a Reference Genome

Alignment maps sequencing reads to positions in a reference genome or transcriptome. The choice of reference and aligner has a major impact on the quality of downstream quantification [6]. RNA-seq aligners must handle reads that span exon-exon junctions, which do not align contiguously to the genome.

Spliced Alignment

STAR is a widely used spliced aligner that can map reads across exon junctions. It is fast and accurate for standard RNA-seq datasets. The aligner requires a genome index built from the reference genome and an annotation file describing gene models. The index construction step is computationally intensive but needs to be performed only once per genome and annotation version.

The choice of reference genome version and annotation is a critical decision. Different versions of the same genome assembly may produce different alignment results, particularly in regions that were updated between versions. The annotation version determines which genes and transcripts are recognized. Analysts should record the exact genome build and annotation version used, because these details are essential for reproducibility and for comparing results across studies.

Mapping Rate Assessment

The mapping rate is the proportion of reads that align to the reference genome. High-quality RNA-seq datasets typically show mapping rates above 70 percent, but the expected rate depends on the organism, the library preparation method, and the quality of the reference genome. Low mapping rates can indicate contamination, poor library quality, or a mismatch between the sample species and the reference genome.

Low mapping rates require investigation before proceeding. Possible causes include adapter contamination that was not removed, ribosomal RNA contamination from incomplete depletion, or sequencing of a different species than expected. Each cause requires a different remedy, so the analyst should diagnose the problem instead of simply lowering alignment stringency.

Multimapping Reads

Some reads align equally well to multiple locations in the genome. These multimapping reads often originate from repetitive regions or gene families with high sequence similarity. The treatment of multimapping reads affects quantification results. Some pipelines discard multimapping reads, while others distribute them among their mapping locations. The choice can influence expression estimates for genes in duplicated or repetitive regions.

Analysts should record the multimapping rate and consider whether it differs systematically between samples or conditions. Consistent multimapping rates across samples are less concerning than variable rates, which may indicate sample-specific technical issues.

Quantification of Gene and Transcript Abundance

Quantification assigns aligned reads to genes or transcripts and produces a count matrix that serves as input for differential expression analysis. The quantification approach determines the resolution of the analysis, with gene-level quantification being the most common and transcript-level quantification providing isoform information.

Gene-Level Quantification with featureCounts

featureCounts is a widely used tool that counts reads overlapping genomic features defined in an annotation file. It is fast, memory-efficient, and suitable for gene-level quantification. The tool requires specification of the feature type, typically exons, and the attribute used to group features into genes, typically the gene identifier.

Read assignment can be ambiguous when a read overlaps multiple genes or when genes overlap on opposite strands. featureCounts offers options for handling multi-overlapping reads, and the analyst must choose a strategy that matches the experimental context. The fraction of reads that are successfully assigned to genes should be recorded for each sample. Low assignment rates indicate problems with annotation, alignment, or library quality.

Strandness and Library Orientation

RNA-seq libraries can be stranded or unstranded, depending on the library preparation protocol. Stranded libraries preserve information about which DNA strand the RNA originated from, which improves the accuracy of quantification for genes with overlapping or antisense transcription. The strandness of the library must be specified correctly for quantification tools, because using the wrong setting will assign reads to the wrong strand and produce incorrect counts.

The strandness can be inferred from the data by examining whether reads align predominantly to the annotated strand of known genes. Many analysis pipelines include a strandness check as part of the quality assessment. Analysts should verify the strandness setting against the library preparation protocol and the observed alignment patterns.

Transcript-Level Quantification

Transcript-level quantification estimates the abundance of individual isoforms, which is more challenging than gene-level quantification because isoforms of the same gene share exonic sequence. Tools such as Salmon and kallisto use lightweight alignment or pseudoalignment approaches to estimate transcript abundances. These methods are fast and can be used for both gene-level and transcript-level analyses.

Isoform-level measurements are less accurate than gene-level measurements, and the accuracy varies by gene and by the degree of isoform diversity [18]. Transcript-level analysis is appropriate when the biological question concerns isoform usage or alternative splicing, but the limitations should be acknowledged in the interpretation of results.

Normalization for Between-Sample Comparisons

Raw read counts are not directly comparable between samples because libraries differ in total sequencing depth and composition. Normalization adjusts counts to account for these technical differences so that expression changes reflect biological variation instead of sequencing artifacts.

Library Size Normalization

The simplest normalization divides counts by the total number of mapped reads in each sample, producing counts per million. This approach assumes that total RNA output is similar across samples, which is not always true. Highly expressed genes in one sample can consume a large fraction of the sequencing capacity and suppress the apparent expression of other genes, a phenomenon known as composition bias.

Median-of-Ratios Normalization in DESeq2

DESeq2 uses a median-of-ratios method that estimates a size factor for each sample based on the ratio of each gene's count to the geometric mean across samples. This approach is robust to composition bias and is appropriate for most differential expression analyses. The size factors are used to scale the counts before statistical testing.

The normalization method should be matched to the downstream statistical analysis. DESeq2, edgeR, and limma-voom each implement normalization approaches that are integrated with their statistical models. Mixing normalization methods and statistical tools without understanding the assumptions can produce misleading results.

Normalization for Other Purposes

Different biological questions require different normalization strategies. For example, comparing expression across genes within a sample requires correction for gene length, because longer genes accumulate more reads. Transcripts per million (TPM) normalization accounts for both gene length and library size and is appropriate for within-sample comparisons. For between-sample comparisons of the same gene, length correction is unnecessary.

Single-cell RNA-seq data require specialized normalization approaches that account for the sparse nature of the data and the technical variability between cells [5][19]. The normalization choices in single-cell analysis differ substantially from bulk RNA-seq and are beyond the scope of this article.

Differential Expression Analysis with DESeq2

Differential expression analysis identifies genes whose expression differs significantly between experimental conditions. DESeq2 is a widely used R package that models count data with a negative binomial distribution and provides statistical tests for comparing groups.

The Statistical Model

DESeq2 estimates the mean and dispersion of counts for each gene and uses a generalized linear model to test for differences between conditions. The model accounts for library size through the size factors estimated during normalization. The dispersion parameter captures the biological variability between replicates, which is essential for accurate hypothesis testing.

The output of DESeq2 includes log2 fold changes, standard errors, test statistics, p-values, and adjusted p-values for each gene. The adjusted p-values control the false discovery rate across the many thousands of genes tested simultaneously. Genes with adjusted p-values below a chosen threshold, commonly 0.05 or 0.1, are considered differentially expressed.

Low-Count Filtering

Genes with very low counts across all samples provide little information and contribute to the multiple testing burden. DESeq2 performs automatic filtering of low-count genes, but analysts can also apply their own filters based on the total counts across samples. The filtering threshold should be chosen to remove genes that cannot be reliably tested while retaining genes with genuine biological signal.

The number of genes passing the low-count filter should be recorded for each analysis. A very small number of detected genes may indicate a problem with library quality, alignment, or quantification, while an unusually large number may suggest contamination or annotation issues.

Interpreting Fold Changes

The log2 fold change reported by DESeq2 represents the estimated magnitude of expression change between conditions. The statistical significance of a gene does not necessarily indicate biological importance, because genes with small fold changes can be statistically significant if their expression is consistent across replicates. Conversely, genes with large fold changes may not reach statistical significance if their expression is highly variable.

Analysts should consider both the magnitude of the fold change and the statistical significance when selecting genes for further study. Many studies apply a fold change threshold in addition to an adjusted p-value threshold, but the choice of thresholds should be justified by the biological context and the goals of the study.

Batch Effects and Confounding

Batch effects are systematic technical differences between groups of samples processed at different times, by different operators, or with different reagent lots. Batch effects can obscure genuine biological differences or create spurious differences if they are confounded with the experimental conditions.

Identifying Batch Effects

Batch effects can be detected by examining the data for patterns that correlate with technical variables instead of biological variables. Principal component analysis is commonly used to visualize the major sources of variation in the data. If the first principal components separate samples by processing batch instead of by experimental condition, a batch effect is likely present.

The experimental design should avoid confounding batch with condition. If all treated samples are processed in one batch and all control samples in another, it is impossible to distinguish treatment effects from batch effects. Randomizing samples across batches is the most effective way to prevent this problem.

Statistical Correction for Known Batches

When batch information is recorded, statistical methods can include batch as a covariate in the differential expression model. DESeq2 supports the inclusion of additional terms in the design formula, which allows the model to account for batch effects while testing for condition effects. This approach is appropriate when the batch assignment is known and the batch effect is consistent across genes.

The design formula should be specified before running the analysis and should reflect the experimental design. Including unnecessary terms reduces statistical power, while omitting important terms can lead to spurious results. The choice of design formula is a key decision that should be documented in the analysis plan.

Handling Unknown or Unmeasured Batch Effects

When batch information is not recorded, or when the batch effect is caused by an unmeasured variable, more complex methods may be needed. Surrogate variable analysis and similar approaches estimate hidden factors from the data and include them in the model. These methods can be useful but require careful application because they can remove genuine biological signal if the biological variation is correlated with the estimated factors.

The best approach to batch effects is prevention through careful experimental design. Analysts should record all relevant technical information about sample processing and include this information in the analysis when possible.

Common Failure Patterns and Troubleshooting

Several problems recur frequently in RNA-seq analysis. Recognizing these patterns and knowing their likely causes can save substantial time and prevent incorrect conclusions.

Low Mapping Rates

A mapping rate well below the expected range for the organism and library type indicates a problem. The first step is to examine the quality control reports for adapter contamination, low-quality bases, and overrepresented sequences. If the reads contain substantial adapter sequence, trimming should be repeated with appropriate parameters. If the reads appear to be high quality but still do not map, the possibility of sample contamination or species mismatch should be investigated.

High Duplication Rates

Excessive PCR duplication can result from over-amplification during library preparation or from low input RNA. Duplicate reads reduce the effective library complexity and can bias expression measurements. If duplication rates are high and variable across samples, the affected samples may need to be re-sequenced or the results interpreted with caution.

Sample Mix-Ups and Contamination

Unexpected expression patterns can indicate that samples were mislabeled or contaminated. Checking the expression of sex-specific genes, species-specific markers, or known housekeeping genes can help identify sample mix-ups. Cross-sample contamination can occur during library preparation and is more difficult to detect. The analysis should include checks for consistency between the observed data and the expected sample characteristics.

Outlier Samples

Individual samples that differ dramatically from others in the same group can distort differential expression results. DESeq2 provides diagnostic plots that help identify outlier samples. The decision to exclude an outlier should be based on evidence of a technical problem instead of on the desire to improve results. Excluding samples without justification reduces the credibility of the analysis.

Discrepancies Between Biological and Technical Replicates

Large variability between biological replicates can reduce statistical power and make it difficult to detect genuine expression differences. The variability may reflect genuine biological heterogeneity or technical problems. Examining the dispersion estimates and the principal component analysis can help distinguish these possibilities.

Reproducibility and Documentation

Reproducibility is a core requirement for RNA-seq analysis. The analysis should be documented in sufficient detail that another researcher can repeat the same steps and obtain the same results. This documentation includes the software versions, parameter settings, reference genome and annotation versions, and the order of analysis steps.

Version Control and Containerization

Analysis pipelines should be managed with version control systems that track changes to scripts and configuration files. Containerization tools such as Docker and Singularity can package software and dependencies in a way that ensures consistent execution across different computing environments. These practices are particularly important for long-term reproducibility, because software updates can change analysis results.

Recording Analysis Parameters

Every parameter choice in the pipeline should be recorded, including trimming thresholds, alignment options, quantification settings, and statistical thresholds. The rationale for each choice should be documented, particularly when the choice deviates from standard practice. This documentation supports both reproducibility and the interpretation of results by other researchers.

Public Data Repositories

Public repositories such as the NCBI provide access to sequencing data and associated metadata [2]. Depositing raw data and analysis results in public repositories supports transparency and enables other researchers to verify or extend the analysis. The NIH Genomic Data Sharing Policy describes expectations for data sharing in NIH-funded research [3]. Researchers should be aware of the data sharing requirements that apply to their funding sources and follow the FAIR principles for making data findable, accessible, interoperable, and reusable [4].

Limitations of RNA-Seq Analysis

RNA-seq provides powerful measurements of gene expression, but the technology and analysis methods have inherent limitations that should be acknowledged in the interpretation of results.

Accuracy of Absolute Measurements

RNA-seq provides accurate measurements of relative expression between genes and between conditions, but it does not provide accurate absolute measurements of transcript abundance [18]. The relationship between read counts and actual transcript numbers is influenced by sequence-specific biases, library preparation efficiency, and other technical factors. Comparisons of expression levels across different genes should be interpreted with this limitation in mind.

Isoform-Level Uncertainty

Quantifying individual transcript isoforms is substantially less accurate than quantifying genes, particularly for genes with many isoforms or isoforms that share large portions of their sequence [18]. Isoform-level results should be validated with orthogonal methods before drawing strong conclusions about isoform switching or alternative splicing.

Detection Limits for Lowly Expressed Genes

Genes with very low expression levels are difficult to measure reliably. The detection limit depends on sequencing depth and library complexity. Genes that are not detected in a sample may be genuinely absent or may be expressed below the detection threshold. The absence of evidence for expression is not evidence for the absence of expression.

Single-Cell Specific Challenges

Single-cell RNA-seq introduces additional analytical challenges, including high data complexity, technical noise, and the need for specialized computational tools [10]. The analysis of single-cell data requires careful quality control at the cell level and the use of methods designed for sparse count data [5][19]. The field lacks universal standardization, which means that different pipelines can produce different results from the same data [10].

Integration with Other Data Types

RNA-seq data are often analyzed in combination with other genomic data types to gain deeper biological insights. The integration of RNA-seq with chromatin accessibility data from ATAC-seq can reveal relationships between chromatin state and gene expression [15]. Multi-omic integration requires careful attention to the different data structures and the appropriate statistical methods for combining them.

Public resources such as GTEx and TCGA provide large collections of RNA-seq data that can be used for secondary analysis and hypothesis generation [12]. These resources have standardized processing pipelines and provide access to processed data, which facilitates comparative analyses. Researchers using public data should understand the processing steps applied to the data and the limitations of the available metadata.

Professional Escalation Criteria

Analysts should recognize when a problem exceeds their expertise and requires consultation with a bioinformatics specialist or a statistician. The following situations warrant escalation:

  • Mapping rates that remain low after troubleshooting standard causes
  • Batch effects that cannot be resolved with available metadata
  • Differential expression results that are highly sensitive to small changes in analysis parameters
  • Discrepancies between RNA-seq results and orthogonal validation experiments
  • Analyses involving non-model organisms with poorly annotated genomes
  • Studies with complex designs involving multiple factors or repeated measures

Early consultation with a specialist is generally more efficient than attempting to resolve unfamiliar problems independently. The cost of incorrect analysis decisions can be substantial, particularly when the results inform further experiments or publications.

Frequently Asked Questions

What is the minimum number of biological replicates needed for differential expression analysis?

The minimum number depends on the biological variability of the system and the magnitude of the expression changes being sought. Statistical methods such as DESeq2 require at least two replicates per condition to estimate biological variability, but three or more replicates per condition are strongly recommended for reliable results. Studies with fewer replicates have low statistical power and may miss genuine differences.

How do I choose between gene-level and transcript-level quantification?

Gene-level quantification is more accurate and is appropriate for most studies that ask which genes change between conditions. Transcript-level quantification is needed when the biological question concerns isoform usage or alternative splicing. Transcript-level measurements are less accurate, and isoform-level results should be interpreted with caution and validated with independent methods.

What should I do if my mapping rate is low?

First examine the quality control reports for adapter contamination, low-quality bases, and overrepresented sequences. Repeat trimming with appropriate parameters if adapters are present. Check whether the reference genome and annotation match the sample species. If the reads are high quality and the reference is correct, consider the possibility of contamination or a library preparation problem.

How can I tell if my data have batch effects?

Perform principal component analysis on the normalized data and examine whether samples cluster by technical variables such as processing date, operator, or reagent lot. If the primary axes of variation separate samples by technical factors instead of biological conditions, a batch effect is likely. Recording detailed metadata about sample processing is essential for detecting and correcting batch effects.

Should I remove duplicate reads from my RNA-seq data?

For standard gene-level differential expression analysis, most pipelines do not remove duplicates by default. Duplicate removal is more important for applications such as variant calling and for libraries with very high duplication rates. The decision depends on the library preparation method and the biological question. Record the duplication rate for each sample and consider whether it varies systematically across conditions.

What is the difference between normalization methods for RNA-seq?

Normalization methods differ in their assumptions about the data. Library size normalization assumes total RNA output is similar across samples. Median-of-ratios normalization used by DESeq2 is robust to composition bias. Transcripts per million normalization accounts for gene length and is appropriate for within-sample comparisons. The choice of normalization should match the statistical analysis and the biological question.

How do I handle genes with zero counts in some samples?

Genes with zero counts in some samples but nonzero counts in others are common in RNA-seq data. DESeq2 and similar tools handle zero counts within their statistical models. Genes with zero counts across all samples provide no information and are typically filtered out before analysis. The filtering threshold should be documented and justified.

Can I compare RNA-seq results across different studies?

Comparing results across studies is challenging because of differences in library preparation, sequencing platforms, analysis pipelines, and annotation versions. Public resources with standardized processing pipelines facilitate cross-study comparisons, but the limitations should be acknowledged. The SEQC project found that measurement performance depends on the platform and analysis pipeline, with large variation for transcript-level profiling [18]. Direct comparisons of absolute expression values across studies are particularly problematic.

Related Bioinformatics Guides

References and Further Reading

This article is educational and does not replace validated analysis plans, institutional policy, clinical interpretation, or specialist review.