Bio applied rna seq analysis

RNA-seq Analysis: Reads to Differential Expression

When to Use

• Planning an RNA-seq experiment (replicate number, sequencing depth, read type, strandedness)
• Deciding between alignment-based (STAR/HISAT2 + featureCounts) and alignment-free (Salmon/kallisto) quantification
• Building or QC-ing a gene x sample count matrix before differential expression
• Converting raw counts to TPM/FPKM for visualization (never for DE testing)
• Running DE testing with DESeq2 (R) or pydeseq2 (Python) and interpreting log2FC/padj

Version Compatibility

STAR ≥2.7.11, HISAT2 ≥2.2.1, Salmon ≥1.10, featureCounts (Subread) ≥2.0, DESeq2 ≥1.42 (Bioconductor ≥3.18), pydeseq2 ≥0.4, Python ≥3.10, R ≥4.3.

Prerequisites

• Python: numpy, pandas, scipy, scikit-learn, pydeseq2; R: DESeq2, tximport (Bioconductor)
• Prior steps: read QC and a GTF/GFF gene annotation for the reference genome

Workflow Overview

FASTQ → QC (FastQC/fastp) → Alignment (STAR/HISAT2) → Counting (featureCounts)
                          → Pseudo-align (Salmon/kallisto) → tximport
                          → Count matrix → Normalization → DESeq2/edgeR → GSEA

Experimental Design Rules

• Replicates: ≥3 biological replicates per condition; more replicates > deeper sequencing
• Depth: 10–30M reads for gene-level DE; 50–100M for isoform discovery
• Read type: SE 50–75 bp for gene-level; PE 100–150 bp for splicing/isoforms
• Strandedness: use strand-specific protocol (dUTP) — necessary for overlapping genes
• Batches: balance conditions across batches; record batch for correction (limma removeBatchEffect or DESeq2 design formula)

Alignment and Quantification

Tool	Memory	Speed	Notes
STAR	30+ GB	Very fast	Most used for human/mouse; outputs GeneCounts
HISAT2	~8 GB	Fast	Lower memory; successor to TopHat2
Salmon/kallisto	<8 GB	10–100x faster	Pseudo-alignment; TPM per transcript

# STAR: build index once, then align + count
STAR --runMode genomeGenerate --genomeDir star_index/ \
     --genomeFastaFiles genome.fa --sjdbGTFfile genes.gtf

STAR --genomeDir star_index/ --readFilesIn R1.fastq R2.fastq \
     --outSAMtype BAM SortedByCoordinate --quantMode GeneCounts

featureCounts -a genes.gtf -o counts.txt -T 4 -p --countReadPairs \
              sample1.bam sample2.bam sample3.bam

# Salmon: alignment-free, ~10-100x faster, comparable gene-level accuracy
salmon index -t transcriptome.fa -i salmon_index
salmon quant -i salmon_index -l A -1 R1.fastq -2 R2.fastq \
             -o sample_quant --validateMappings

Salmon/kallisto output TPM per transcript — aggregate to gene level with tximport (R) before DESeq2, or sum manually in Python.

Normalization

Unit	Formula	Cross-sample comparable?	Use for DE?
Raw counts	—	No	Yes (DESeq2/edgeR)
RPKM/FPKM	counts / (len_kb × total_M)	No	No
TPM	(counts/len) / sum(counts/len) × 1e6	Yes (sums to 1M)	No

Goal: compute TPM for reporting/visualization and DESeq2-style size factors for DE-ready normalization. Approach: length-normalize then scale to a million for TPM; use the median-of-ratios method (geometric mean per gene, then per-sample median ratio) for size factors, since it is robust to a few highly-expressed genes.

import numpy as np
import pandas as pd


def counts_to_tpm(counts, lengths_bp):
    """Convert a 1D array of raw counts to TPM. lengths_bp is gene length in bp."""
    rate = counts / lengths_bp
    return rate / rate.sum() * 1e6


def deseq2_size_factors(count_matrix: pd.DataFrame) -> pd.Series:
    """Compute DESeq2-style size factors (median-of-ratios method).

    count_matrix: genes (rows) x samples (columns), raw integer counts.
    """
    nonzero_mask = (count_matrix > 0).all(axis=1)
    filtered = count_matrix.loc[nonzero_mask]

    log_means = np.log(filtered).mean(axis=1)
    geo_means = np.exp(log_means)

    size_factors = {}
    for sample in filtered.columns:
        ratios = filtered[sample] / geo_means
        size_factors[sample] = np.median(ratios)
    return pd.Series(size_factors)

Rule: use raw counts for DE testing; use TPM only for visualization/reporting — TPM/FPKM already normalize in a way that breaks count-based negative-binomial models.

Count Matrix

• Shape: (n_genes × n_samples), integer values, negative-binomial (overdispersed Poisson) distributed
• Zeros are common for lowly expressed genes — do not impute
• Filter low-count genes before DE (e.g., keep genes with ≥10 counts in ≥3 samples) to reduce multiple-testing burden

import pandas as pd

# Load featureCounts output; drop its 5 annotation meta-columns, keep sample columns
counts = pd.read_csv('counts.txt', sep='\t', comment='#', index_col=0)
counts = counts.iloc[:, 5:]

Differential Expression

Goal: test each gene for a mean-expression difference between conditions, controlling the false discovery rate. Approach: DESeq2 fits a negative-binomial GLM per gene with shrunk dispersion estimates and a Wald test, then applies Benjamini-Hochberg correction. pydeseq2 reimplements the same statistical model in Python.

library(DESeq2)

dds <- DESeqDataSetFromMatrix(
    countData = count_matrix,
    colData   = sample_info,
    design    = ~ condition
)
dds <- DESeq(dds)  # size factors, dispersion estimation, Wald test

res <- results(dds, contrast = c("condition", "Treatment", "Control"))
res_shrunk <- lfcShrink(dds, coef = "condition_Treatment_vs_Control", type = "apeglm")
summary(res)

from pydeseq2.dds import DeseqDataSet
from pydeseq2.ds import DeseqStats

dds = DeseqDataSet(counts=counts_df, metadata=sample_info, design_factors="condition")
dds.deseq2()

stat_res = DeseqStats(dds, contrast=["condition", "Treatment", "Control"])
stat_res.summary()
res = stat_res.results_df  # log2FoldChange, pvalue, padj

Pitfalls

• Coordinate systems: BED is 0-based; GTF/GFF are 1-based — off-by-one errors when mixing
• Batch effects: PCA or hierarchical clustering before DE to detect; correct with design formula or removeBatchEffect — never correct counts and feed to DESeq2
• Multiple testing: always use padj (Benjamini-Hochberg); never filter on raw p-value alone
• Library size confounding: do not use TPM/FPKM as input to DESeq2/edgeR — they have already normalized in a way that breaks count-based models
• Poly-A vs rRNA depletion: poly-A misses non-coding RNA and degraded RNA; use rRNA depletion for FFPE or total RNA experiments
• Strandedness mismatch: wrong --strandedness in featureCounts halves your counts; always verify with RSeQC infer_experiment.py
• Outlier samples: one bad sample can dominate DE results; always run PCA/hclust QC before DE

See Also

• rnaseq — DESeq2/edgeR normalization, PCA, volcano plots, GSEA in more depth
• bio-applied-mirna-seq-pipeline — small RNA-seq variant of this workflow
• bio-applied-scrna-preprocessing — single-cell analog (no bulk count matrix)
• bio-applied-isoform-analysis — transcript-level/splicing follow-up analysis