RNA-seq Data Processing: A Step-by-Step Guide

By Dr. Zubair Khalid, DVM, MS, PhD · 2026-08-11

Introduction to RNA-seq Data Processing

What is RNA-seq data processing?

RNA sequencing (RNA-seq) data processing is the computational pipeline that converts raw sequencing reads—the output of a high-throughput sequencing instrument—into interpretable biological measurements of gene and transcript expression. The raw data from an Illumina sequencer consists of millions to billions of short nucleotide sequences (typically 50–150 base pairs) stored in FASTQ format, each accompanied by a per-base quality score. These reads are meaningless without processing: they must be assessed for quality, cleaned of technical artifacts, aligned to a reference, and quantified to produce a count matrix where rows are genes or transcripts and columns are biological samples.

The purpose of RNA-seq data processing is to eliminate technical variation introduced during library preparation and sequencing while preserving biological signal. This distinction matters because technical noise—adapter contamination, sequencing errors, PCR duplicates, and mapping artifacts—can obscure or mimic differential expression if left unaddressed. A well-processed dataset yields accurate expression estimates that serve as the foundation for differential expression analysis, isoform discovery, and functional interpretation.

Overview of the RNA-seq workflow

The standard RNA-seq processing workflow follows a logical sequence of steps, each building on the previous:

1. Quality control (QC) of raw FASTQ files to assess read quality, GC content, adapter contamination, and duplication levels.
2. Preprocessing to trim low-quality bases and adapter sequences, and optionally filter reads entirely.
3. Alignment of cleaned reads to a reference genome or transcriptome using a splice-aware aligner.
4. Quantification of expression at the gene or transcript level using count-based tools or pseudo-alignment algorithms.
5. Normalization to make samples comparable by accounting for sequencing depth and library composition.
6. Differential expression analysis using statistical models that estimate variance and test for significant changes.
7. Downstream analysis including dimensionality reduction, clustering, and functional enrichment.

Each step has multiple tool options, and the choices made at each stage propagate through the pipeline. The sections that follow provide mechanistic detail on each step, including the algorithms underlying common tools and the practical decisions you will face.

Quality Control and Preprocessing of Raw Reads

Assessing read quality with FastQC

The first step in any RNA-seq processing pipeline is to assess the quality of your raw sequencing data. FastQC is the de facto standard for this purpose. It generates a comprehensive HTML report containing 12 modules that evaluate different aspects of read quality, including per-base quality scores, per-sequence quality scores, GC content distribution, adapter content, and overrepresented sequences.

The per-base quality plot shows the Phred quality score (Q score) for each position across all reads. A Phred score of Q30 corresponds to a base call accuracy of 99.9% (one error per 1,000 bases), while Q20 corresponds to 99% accuracy. For Illumina sequencing, you should expect median quality scores above Q30 for most positions. A sharp drop in quality toward the 3′ end of reads is normal, particularly for longer read lengths, but a precipitous decline in the first 10–20 bases may indicate problems with the sequencing run or library preparation.

The GC content distribution should approximate a normal distribution centered near the genomic GC content of your organism (approximately 41% for human). A bimodal distribution may indicate contamination from another species, while a sharp spike at a specific GC percentage often signals adapter dimers or PCR artifacts. Overrepresented sequences—those appearing more than 0.1% of the time—frequently represent adapter sequences, rRNA contamination, or highly expressed transcripts.

FastQC is a diagnostic tool, not a filtering tool. It tells you what problems exist but does not fix them. The decision of whether to trim, filter, or re-sequence depends on the severity of the issues identified. As a rule of thumb, if median per-base quality drops below Q20 at any position, or if adapter content exceeds 5% of reads, trimming is warranted.

Trimming adapters and low-quality bases

Adapter trimming is essential for RNA-seq because fragment lengths are often shorter than the read length. When a fragment is shorter than the sequencing read, the sequencer continues reading into the adapter sequence, producing reads that contain adapter sequence at the 3′ end. If these adapters are not removed, they will cause alignment failures or mismapped reads.

Two widely used tools for trimming are Trimmomatic and cutadapt. Both operate on the same principle: identify adapter sequences in reads and remove them, along with low-quality bases.

Trimmomatic uses a sliding-window approach. The command:

trimmomatic PE -phred33 input_R1.fastq.gz input_R2.fastq.gz \
  output_R1_paired.fastq.gz output_R1_unpaired.fastq.gz \
  output_R2_paired.fastq.gz output_R2_unpaired.fastq.gz \
  ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 \
  LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36

performs four operations: (1) ILLUMINACLIP removes adapter sequences allowing up to 2 mismatches, requiring a 30-base pair match for read-through and a 10-base pair match for palindrome mode in paired-end reads; (2) LEADING:3 and TRAILING:3 remove bases with quality below Q3 at the 5′ and 3′ ends, respectively; (3) SLIDINGWINDOW:4:15 scans a 4-base window and trims once the average quality drops below Q15; and (4) MINLEN:36 discards reads shorter than 36 bases after trimming.

cutadapt operates similarly but uses a more flexible adapter-matching algorithm based on semi-global alignment. It can handle 3′ adapters, 5′ adapters, and linked adapters, making it more versatile for complex library preparations. A typical command:

cutadapt -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCA \
  -q 20 -m 36 -o output.fastq.gz input.fastq.gz

removes the Illumina TruSeq adapter sequence, trims bases with quality below Q20, and discards reads shorter than 36 bases.

The choice between Trimmomatic and cutadapt is largely personal preference; both produce comparable results. The key parameters to adjust are the quality threshold (Q15–Q20 is typical), the minimum read length (30–36 bases is common), and the adapter sequences, which depend on your library preparation kit.

Post-trimming QC checks

After trimming, you should rerun FastQC on the cleaned reads to confirm that adapter contamination has been removed and quality profiles have improved. You should also check the proportion of reads retained after trimming—losing more than 20–30% of reads suggests either poor initial data quality or overly aggressive trimming parameters.

An additional QC metric at this stage is the read duplication rate. High duplication levels (above 50%) can indicate PCR over-amplification during library preparation, which biases expression estimates. However, in RNA-seq, high duplication rates are expected for highly expressed genes, so this metric should be interpreted with caution. Tools like Picard MarkDuplicates can identify and flag PCR duplicates, but for standard RNA-seq analysis, duplicate removal is generally not recommended because it can remove legitimate biological reads from highly expressed genes.

Read Alignment to a Reference Genome

Choosing an aligner: STAR vs. HISAT2

Alignment is the process of determining where each read originated in the genome. For RNA-seq, this requires a splice-aware aligner—one that can split reads across exon-exon junctions. Two aligners dominate the field: STAR (Spliced Transcripts Alignment to a Reference) and HISAT2 (Hierarchical Indexing for Spliced Alignment of Transcripts 2).

STAR uses a two-step approach based on a maximal mappable prefix (MMP) algorithm. In the first step, STAR finds the longest sequence that matches the genome exactly, allowing for mismatches. When it encounters a mismatch or the end of the read, it extends the search from the next base. This allows STAR to find the MMP for each read segment. In the second step, STAR stitches together the MMPs to form the full read alignment, allowing for gaps that represent introns. STAR's seed-and-extend approach is extremely fast—it can align 100 million paired-end reads in under an hour on a standard server—but it requires substantial memory (approximately 30 GB for the human genome).

HISAT2 uses a different strategy based on a hierarchical graph FM index (HGFM). It builds a global FM index of the entire genome and a set of local FM indexes for smaller genomic regions. This hierarchical approach allows HISAT2 to achieve speed comparable to STAR while using significantly less memory (approximately 5 GB for the human genome). HISAT2 also supports alignment to a graph genome that includes known genetic variants, which can improve alignment accuracy for samples with common polymorphisms.

For most applications, both aligners produce comparable results. STAR is generally preferred for its accuracy and speed, particularly for large datasets, while HISAT2 is a better choice when memory is constrained. A typical STAR alignment command:

STAR --runMode alignReads \
  --genomeDir /path/to/STAR_index \
  --readFilesIn sample_R1.fastq.gz sample_R2.fastq.gz \
  --readFilesCommand zcat \
  --outSAMtype BAM SortedByCoordinate \
  --outFileNamePrefix sample_ \
  --quantMode GeneCounts

The --quantMode GeneCounts option generates per-gene counts directly from the alignment, which can be used for downstream analysis without a separate quantification step.

Handling spliced reads

The critical feature of RNA-seq aligners is their ability to handle reads that span exon-exon junctions. A read of 100 bases that spans a junction might have 60 bases mapping to the end of one exon and 40 bases mapping to the beginning of the next exon, with the intervening intronic sequence absent from the read.

STAR handles this through its MMP algorithm: it finds the maximal match in the first exon, then continues searching from the next base in the second exon, ultimately reporting the alignment with a gap (the intron) between the two segments. HISAT2 uses its graph FM index to find spliced alignments by searching for the read sequence across exon boundaries.

Both aligners use annotated splice junctions to guide alignment. The genome index includes information about known exon-exon junctions from the annotation file (GTF/GFF format), which improves alignment accuracy for reads spanning known junctions. However, both aligners can also discover novel junctions—splice sites not present in the annotation—which is essential for detecting new isoforms.

The output of alignment is a SAM (Sequence Alignment/Map) file, or its compressed binary form, BAM. Each line in a SAM file represents one read alignment and contains 11 mandatory fields, including the read name, flag (which encodes strand and pairing information), reference name, position, mapping quality (MAPQ), and the CIGAR string that describes the alignment (e.g., 60M200N40M indicates 60 bases matched, 200 bases skipped (intron), 40 bases matched).

Generating and filtering alignments (SAM/BAM files)

After alignment, you should filter the BAM file to remove poor-quality alignments. Common filters include:

• Mapping quality (MAPQ): Remove reads with MAPQ below 10–20. MAPQ is a Phred-scaled probability that the read is incorrectly mapped; a MAPQ of 20 means a 1% chance of misalignment.
• Secondary alignments: Reads that map to multiple locations (multi-mappers) receive secondary alignment flags. For most analyses, you should retain only primary alignments (flag 0 or 16) and discard secondary (flag 256) and supplementary (flag 2048) alignments.
• Improper pairs: For paired-end data, remove read pairs that are not properly paired (flag 2 not set), as these may represent chimeric fragments or alignment errors.

The samtools suite provides the standard tools for BAM file manipulation:

samtools view -b -q 20 -F 0x100 -F 0x800 sample_sorted.bam > sample_filtered.bam
samtools sort -o sample_sorted.bam sample_filtered.bam
samtools index sample_sorted.bam

The -q 20 flag removes reads with MAPQ below 20, and the -F flags remove secondary and supplementary alignments. After filtering, you should check the alignment rate—the proportion of reads that mapped to the genome. For human RNA-seq data, alignment rates of 80–90% are typical; lower rates may indicate contamination, poor reference genome quality, or issues with library preparation.

Quantifying Gene and Transcript Expression

Gene-level quantification with featureCounts

Gene-level quantification counts the number of reads that overlap each gene in the annotation. featureCounts (from the Subread package) is the most widely used tool for this purpose. It assigns reads to genes based on the genomic coordinates of the gene's exons, using the annotation file (GTF format) to define gene models.

featureCounts uses a chromosomal sweeping algorithm that processes reads in genomic order, which makes it highly efficient. It counts reads that overlap exons of each gene, with several options for handling reads that overlap multiple genes or features:

• -s (strandness): Specify whether the library is unstranded (0), stranded (1), or reversely stranded (2). This is critical for accurate counting—using the wrong strand setting can result in counting reads from the wrong strand and severely biasing expression estimates.
• -p: Specify that the data is paired-end. For paired-end data, featureCounts counts fragments (pairs of reads) rather than individual reads, which avoids double-counting.
• -t and -g: Specify the feature type (default: exon) and attribute type (default: gene_id) to use for grouping.

A typical command:

featureCounts -a annotation.gtf -o counts.txt \
  -s 2 -p --countReadPairs \
  sample1_sorted.bam sample2_sorted.bam sample3_sorted.bam

The output is a count matrix where each row is a gene and each column is a sample, with the number of reads (or fragments) assigned to each gene.

HTSeq-count is an alternative gene-level counter that uses a more conservative approach: it only counts reads that map uniquely to a single gene and discards reads that overlap multiple genes (ambiguous reads). This makes HTSeq-count more stringent but also more likely to underestimate expression of genes with overlapping or closely spaced annotations. featureCounts is generally preferred for its speed and flexibility.

Transcript-level quantification with Salmon/kallisto

Gene-level quantification is sufficient for most differential expression analyses, but it cannot distinguish between isoforms of the same gene. Transcript-level quantification estimates the abundance of each transcript isoform, which is essential for studying alternative splicing or isoform switching.

Two tools dominate transcript-level quantification: Salmon and kallisto. Both use pseudo-alignment or lightweight alignment approaches that bypass full genomic alignment. Instead of finding the exact genomic position of each read, these tools determine which transcripts the read is compatible with, based on the sequence of k-mers (short subsequences of length k, typically 31) in the read.

kallisto uses a technique called pseudo-alignment based on a colored de Bruijn graph. It builds a graph where nodes represent k-mers from the transcriptome, and edges connect k-mers that are adjacent in transcripts. Reads are then "pseudo-aligned" by finding the path through the graph that matches the read sequence, which identifies the set of transcripts compatible with the read. This approach is extremely fast—kallisto can process 10 million reads in under 5 minutes on a standard laptop.

Salmon uses a similar k-mer-based approach but with a more sophisticated model that accounts for sequence-specific biases (e.g., GC bias, position-specific bias) and fragment length distributions. Salmon also supports selective alignment, which uses a standard aligner (like RapMap) to find candidate alignment positions before quantifying, providing a middle ground between full alignment and pseudo-alignment.

Both tools use an expectation-maximization (EM) algorithm to estimate transcript abundances. The EM algorithm iteratively assigns reads to transcripts based on the current abundance estimates, then updates the estimates based on these assignments, converging to maximum-likelihood estimates of transcript abundance. The output is expressed in transcripts per million (TPM), which normalizes for both sequencing depth and transcript length.

A typical Salmon command:

salmon quant -i salmon_index -l A \
  -1 sample_R1.fastq.gz -2 sample_R2.fastq.gz \
  -p 8 --validateMappings -o sample_quant

The -l A flag tells Salmon to automatically determine the library type (strandedness), which is useful when you are unsure of the library preparation protocol.

Understanding count matrices

Regardless of the quantification tool, the output is a count matrix. For gene-level quantification, each entry is an integer count representing the number of reads or fragments assigned to a gene in a sample. For transcript-level quantification, each entry is an estimated count (which may be fractional due to the EM algorithm) or TPM.

Count matrices are the input for downstream statistical analysis. The key properties of count data are:

• Discreteness: Counts are integers (or nearly so), which means they follow a negative binomial distribution rather than a normal distribution.
• Heteroscedasticity: The variance of counts increases with the mean—highly expressed genes have higher variance than lowly expressed genes.
• Dynamic range: Counts can range from 0 to hundreds of thousands, spanning several orders of magnitude.

These properties dictate the choice of statistical methods for normalization and differential expression analysis, as discussed in the following sections.

Normalization Methods for Between-Sample Comparisons

Why normalization is necessary

Raw read counts are not directly comparable between samples because of differences in sequencing depth (total number of reads per sample) and library composition (the relative proportions of transcripts). A sample sequenced to 50 million reads will have roughly twice the counts of a sample sequenced to 25 million reads for the same gene, even if expression is identical. Similarly, if one sample has a small number of highly expressed genes that consume a large fraction of the reads, all other genes will appear less expressed in that sample, even if their absolute expression is unchanged.

Normalization methods adjust raw counts to account for these technical factors, making expression values comparable between samples.

Common normalization methods

Several normalization methods are in common use, each with different assumptions and applications:

RPKM and FPKM are calculated as:

RPKM = (reads mapped to gene × 10^9) / (total mapped reads × gene length in kb)

TPM is calculated as:

TPM = (reads per transcript length) / (sum of all reads per transcript length) × 10^6

The key difference between RPKM/FPKM and TPM is the order of normalization. RPKM normalizes for length first, then for sequencing depth, while TPM normalizes for length first, then scales to a constant sum (1 million) across all transcripts. This makes TPM more consistent across samples and is the preferred metric for comparing expression levels within and between samples.

Median-of-ratios (used by DESeq2) and TMM (used by edgeR) are more sophisticated methods designed specifically for differential expression analysis. Instead of normalizing by total reads, they estimate size factors that account for library composition differences. The median-of-ratios method calculates the geometric mean of counts for each gene across all samples, then computes the ratio of each sample's count to this geometric mean, and uses the median of these ratios as the size factor. TMM uses a similar approach but trims the most extreme ratios and uses a weighted average.

When to use which method

The choice of normalization method depends on your analysis goal:

• For differential expression analysis with DESeq2 or edgeR, use the built-in normalization (median-of-ratios or TMM, respectively). These methods are designed to minimize false positives caused by composition biases.
• For comparing expression levels within a sample (e.g., which genes are most highly expressed), use TPM.
• For cross-sample comparisons of absolute expression (e.g., comparing expression of a gene across tissues), TPM is appropriate, but be aware that it does not correct for composition biases.
• Avoid RPKM/FPKM for cross-sample comparisons, as they are mathematically inconsistent between samples. The sum of RPKM values across genes differs between samples, which makes them unreliable for comparing expression levels across samples.

Differential Expression Analysis

Designing the model matrix

Differential expression analysis identifies genes whose expression differs significantly between experimental conditions. The standard tools—DESeq2, edgeR, and limma-voom—all use generalized linear models (GLMs) to model read counts as a function of experimental design.

The first step is to define the model matrix, which encodes the experimental design. For a simple two-group comparison (e.g., treated vs. control), the model matrix has two columns: an intercept (representing the control group) and a coefficient for the treatment effect. For more complex designs (e.g., time course, multiple factors, paired samples), the model matrix includes additional columns for each factor and interaction term.

In DESeq2, the design formula is specified as:

dds <- DESeqDataSetFromMatrix(countData = counts, 
                              colData = metadata, 
                              design = ~ condition)

The design formula ~ condition tells DESeq2 to model counts as a function of the condition factor. For a paired design (e.g., same patient before and after treatment), you would use ~ patient + condition to account for patient-to-patient variation.

Dispersion estimation and shrinkage

A critical challenge in differential expression analysis is that RNA-seq experiments typically have few biological replicates (often 3–5 per condition), which makes it difficult to estimate the variance of each gene's expression accurately. The negative binomial distribution models count data with two parameters: the mean and the dispersion (which quantifies the variance relative to the mean).

DESeq2 and edgeR use a two-step approach to estimate dispersion:

1. Per-gene dispersion estimation: For each gene, the dispersion is estimated from the observed counts using maximum likelihood.
2. Dispersion shrinkage: The per-gene estimates are shrunk toward a fitted curve that models dispersion as a function of mean expression. This empirical Bayes approach borrows information from genes with similar expression levels to stabilize estimates for genes with high variance or low counts.

The shrinkage is particularly important for genes with low read counts, where the dispersion estimate is noisy. Without shrinkage, these genes would have inflated variance estimates, reducing statistical power to detect differential expression.

In DESeq2, dispersion shrinkage is automatic. The DESeq() function performs the full pipeline:

dds <- DESeq(dds)
res <- results(dds)

The results() function extracts the results table containing the log2 fold change, standard error, test statistic, and adjusted p-value for each gene.

Multiple testing correction (FDR)

When testing thousands of genes simultaneously, the probability of false positives increases dramatically. If you test 20,000 genes at a significance threshold of p < 0.05, you would expect 1,000 false positives by chance alone. Multiple testing correction controls the number of false positives.

The most common approach is the Benjamini-Hochberg (BH) procedure, which controls the false discovery rate (FDR)—the expected proportion of false positives among all rejected hypotheses. The BH procedure ranks genes by p-value, then applies a threshold that depends on the rank and total number of tests. The result is an adjusted p-value (padj) for each gene.

In DESeq2, the default multiple testing correction is the BH procedure, and the results table includes the padj column. A common threshold for significance is padj < 0.05, which means that approximately 5% of the genes identified as differentially expressed are expected to be false positives.

The number of biological replicates has a direct impact on statistical power. With 3 replicates per condition, you can reliably detect large fold changes (≥2-fold) for moderately expressed genes. With 5–10 replicates, you can detect smaller fold changes and identify differential expression for lowly expressed genes. Increasing replication is almost always more beneficial than increasing sequencing depth once you have at least 10–20 million reads per sample.

Downstream Analysis and Visualization

Sample-level QC with PCA

Before interpreting differential expression results, you should perform sample-level quality control to confirm that biological replicates cluster together and that samples separate by experimental condition. Principal component analysis (PCA) is the standard tool for this purpose.

PCA reduces the high-dimensional expression data (thousands of genes) to a few principal components that capture the largest sources of variation. The first principal component (PC1) captures the direction of maximum variance, PC2 captures the second-largest variance orthogonal to PC1, and so on. Plotting PC1 versus PC2 reveals the overall structure of the data: samples with similar expression profiles cluster together, while samples with distinct profiles separate.

In DESeq2, you can generate a PCA plot using the variance-stabilizing transformation (VST) or the regularized log (rlog) transformation, which stabilize the variance across the range of expression values:

__MASK_11__

A well-behaved dataset shows clear separation between conditions and tight clustering within conditions. If replicates do not cluster, this may indicate batch effects, sample mislabeling, or technical problems that need to be addressed before proceeding with downstream analysis.

Visualizing expression with heatmaps

Heatmaps provide a visual summary of expression patterns across samples and genes. The pheatmap R package is commonly used for this purpose. A typical heatmap displays genes as rows and samples as columns, with colors representing expression levels (usually z-scored across samples).

For visualizing differentially expressed genes, you typically select the top genes by significance (e.g., padj < 0.05 and |log2 fold change| > 1) and display their expression across all samples. Heatmaps are particularly useful for identifying co-expressed gene clusters and for confirming that the expression patterns align with the experimental design.

Functional enrichment analysis

Once you have a list of differentially expressed genes, the next question is: what biological processes or pathways are these genes involved in? Gene ontology (GO) enrichment analysis tests whether your gene list is enriched for specific functional categories compared to the background set of all genes.

The standard tool for GO enrichment is clusterProfiler in R, which uses a hypergeometric test to identify overrepresented GO terms:

library(clusterProfiler)
library(org.Hs.eg.db)

ego <- enrichGO(gene = significant_genes,
                OrgDb = org.Hs.eg.db,
                ont = "BP",
                pAdjustMethod = "BH",
                qvalueCutoff = 0.05)

The output lists GO terms ranked by enrichment significance, along with the gene ratio (proportion of your gene list annotated to the term) and the adjusted p-value.

KEGG pathway enrichment is a complementary approach that tests for enrichment of specific signaling or metabolic pathways. The same clusterProfiler package provides enrichKEGG() for this purpose.

A common pitfall in enrichment analysis is the selection bias introduced by choosing a significance threshold. Genes with high expression are more likely to be detected as differentially expressed, so your gene list may be biased toward highly expressed genes, which in turn may be enriched for certain functional categories. Using a less stringent threshold (e.g., padj < 0.1) or ranking-based approaches (e.g., GSEA) can mitigate this bias.

Common Pitfalls and Best Practices in RNA-seq Processing

Strandness and library prep

One of the most common sources of error in RNA-seq processing is specifying the wrong strandness (also called library type). RNA-seq libraries can be:

• Unstranded: Reads map to both strands of the genome, and the strand of the read does not indicate the strand of the original RNA molecule.
• Stranded (forward): The read maps to the same strand as the RNA (i.e., the sense strand).
• Reversely stranded (reverse): The read maps to the opposite strand of the RNA (i.e., the antisense strand).

The strandness of your library depends on the library preparation kit. For example, the Illumina TruSeq Stranded mRNA kit produces reversely stranded libraries, while the NEBNext Ultra Directional RNA Library Prep Kit also produces reversely stranded libraries. If you specify the wrong strandness in featureCounts (-s parameter) or in Salmon (-l parameter), you will count reads from the wrong strand, leading to severe underestimation of expression for genes on the opposite strand.

To determine strandness empirically, you can use tools like infer_experiment.py from the RSeQC package, which examines the strand of reads relative to annotated genes in a small subset of your aligned data.

Reference genome version and annotation

The choice of reference genome version and annotation file has a major impact on your results. Different versions of the human genome (e.g., GRCh37 vs. GRCh38) differ in sequence content and gene annotations, and using mismatched versions can lead to alignment errors and incorrect gene counts.

Best practices:

• Use the same genome version and annotation for all samples in a study.
• Download the annotation file (GTF) that matches the genome version exactly. The Ensembl and GENCODE annotations are updated regularly, and using an outdated GTF with a newer genome can cause errors.
• For human data, GENCODE annotation is generally preferred for its comprehensive gene models.
• Document the genome version and annotation file version in your methods section, as this is essential for reproducibility.

Reproducibility and documentation

RNA-seq processing involves many steps, each with multiple parameters. To ensure reproducibility:

• Use a workflow management system such as Snakemake, Nextflow, or CWL to define the entire pipeline in a single file. This makes it easy to rerun the analysis and to share the pipeline with collaborators.
• Record software versions for all tools used. Version differences can affect results, particularly for aligners and quantification tools.
• Store parameter settings in configuration files rather than hard-coding them in commands.
• Save intermediate files (trimmed reads, BAM files, count matrices) so that you can trace any issues back to the specific step where they occurred.

The Process Validation principles used in manufacturing—defining the process, qualifying the inputs, and verifying the outputs—apply equally well to computational pipelines. Similarly, the documentation standards used in FDA Approval Process for Biologics provide a model for the level of detail required for reproducible computational analysis.

Frequently Asked Questions

How do I process RNA-seq data from raw FASTQ files?

The standard workflow is: (1) run FastQC on raw FASTQ files to assess quality; (2) trim adapters and low-quality bases with Trimmomatic or cutadapt; (3) rerun FastQC to confirm quality improvement; (4) align reads to a reference genome with STAR or HISAT2; (5) quantify gene expression with featureCounts or transcript expression with Salmon/kallisto; (6) normalize counts; and (7) perform differential expression analysis with DESeq2 or edgeR. Each step is described in detail in the sections above.

What is the best RNA-seq data processing pipeline?

There is no single "best" pipeline—the optimal choice depends on your experimental design, computational resources, and analysis goals. For most applications, a pipeline using STAR for alignment and featureCounts for gene-level quantification, followed by DESeq2 for differential expression, is robust and well-documented. For transcript-level analysis or when computational resources are limited, Salmon or kallisto provide faster alternatives with comparable accuracy.

How do I choose between alignment-based and pseudo-alignment tools?

Alignment-based tools (STAR, HISAT2) provide a BAM file that can be used for additional analyses such as variant calling, splice junction discovery, or visualization in a genome browser. Pseudo-alignment tools (Salmon, kallisto) are faster and require less memory, but they do not produce alignments suitable for these downstream applications. If you only need expression quantification, pseudo-alignment is sufficient; if you need alignments for other purposes, use an alignment-based approach.

What is the difference between RPKM, FPKM, and TPM?

RPKM (reads per kilobase per million) and FPKM (fragments per kilobase per million) normalize read counts by gene length and sequencing depth. FPKM is the paired-end equivalent of RPKM. TPM (transcripts per million) normalizes by transcript length first, then scales to a constant sum of 1 million across all transcripts. TPM is mathematically more consistent for cross-sample comparisons and is preferred over RPKM/FPKM.

How do I know if my RNA-seq data is stranded?

You can determine strandness empirically using the infer_experiment.py script from RSeQC, which compares the strand of aligned reads to annotated genes. Alternatively, you can check the documentation for your library preparation kit—most commercial kits specify whether they produce stranded or unstranded libraries.

What are common mistakes in RNA-seq data processing?

Common mistakes include: specifying the wrong strandness, using mismatched genome and annotation versions, skipping quality control, over-trimming reads, using too few biological replicates, and applying the wrong normalization method for the analysis goal.

How long does it take to process RNA-seq data?

Processing time depends on the number of samples, read depth, and computational resources. For a typical dataset of 10 samples with 30 million paired-end reads each, alignment with STAR takes approximately 1–2 hours on a server with 16 cores, quantification with featureCounts takes minutes, and differential expression analysis with DESeq2 takes minutes. Pseudo-alignment with Salmon is substantially faster, processing all 10 samples in under 30 minutes.

Key Takeaways

• RNA-seq data processing transforms raw FASTQ reads into interpretable expression measurements through a pipeline of quality control, trimming, alignment, quantification, normalization, and statistical analysis.
• Quality control with FastQC before and after trimming is essential for identifying adapter contamination, low-quality bases, and other technical artifacts that can bias results.
• Splice-aware aligners (STAR, HISAT2) are required for RNA-seq because reads span exon-exon junctions; the choice of aligner depends on computational resources and downstream analysis needs.
• Gene-level quantification with featureCounts is sufficient for most differential expression analyses, while transcript-level quantification with Salmon or kallisto is required for isoform-level studies.
• Normalization is critical for between-sample comparisons; TPM is appropriate for comparing expression levels, while DESeq2's median-of-ratios and edgeR's TMM are designed for differential expression analysis.
• Differential expression analysis with DESeq2 or edgeR uses negative binomial models with dispersion shrinkage and multiple testing correction to identify statistically significant changes.
• Always document software versions, parameters, and reference genome versions, and use workflow management systems to ensure reproducibility of your analysis.

Related Topics

• RNA-seq vs Scrna seq
• Read ChIP-seq Data
• Single Nuclear RNA-seq
• Sanger Sequencing Protocol
• Bisulfite Sequencing

Related Clinical & Scientific Guides

• MAPK Pathway: Mechanism, Function, and Clinical Relevance
• Mammalian Cell Culture Bioreactors: A Practical Guide
• Nucleotide Formation: Biosynthesis and Assembly of DNA/RNA Building Blocks