rnaseq-analysis

name rnaseq-analysis description Run bulk RNA-seq differential expression from a gene x sample count matrix using DESeq2 (R), pydeseq2, or edgeR — normalization (median-of-ratios/TPM), Wald/LRT testing, BH-adjusted p-values, volcano/MA plots, and GSEA/ORA. Use when doing RNA-seq DE, comparing treatment vs control expression, building a FASTQ-to-DESeq2 pipeline, or asked about TPM/RPKM/FPKM, count matrices, log2FoldChange, padj, or STAR/Salmon/featureCounts/tximport. tool_type r primary_tool DESeq2 RNA-seq Differential Expression Analysis When to Use Testing which genes differ between conditions (treatment vs. control, tumor vs. normal) from RNA-seq count data Deciding between alignment-based (STAR + featureCounts) and alignment-free (Salmon/kallisto + tximport) quantification Choosing/justifying a normalization unit (raw counts vs. TPM/RPKM/FPKM) for a given use case Building volcano plots, MA plots, or PCA/clustering QC from a count matrix Following up DE results with gene set enrichment (GSEA/ORA) Version Compatibility DESeq2 ≥1.42 (Bioconductor 3.18+), apeglm ≥1.24, edgeR ≥4.0, R ≥4.3, pydeseq2 ≥0.4, Python ≥3.10, pandas ≥2.0, numpy ≥1.26, scipy ≥1.11, scikit-learn ≥1.3, statsmodels ≥0.14. Prerequisites R: BiocManager::install(c("DESeq2","apeglm")) . Python: pip install pydeseq2 pandas numpy scipy scikit-learn statsmodels matplotlib seaborn . CLI tools (if starting from FASTQ): STAR or HISAT2, Subread (featureCounts), Salmon, and tximport (R) for alignment-free counts. Concepts: count matrix layout (genes x samples), negative-binomial count model, Benjamini-Hochberg FDR. Workflow FASTQ -> QC (FastQC/fastp) -> Alignment (STAR/HISAT2) -> featureCounts/HTSeq OR -> Pseudoalignment (Salmon/kallisto) -> tximport -> Count matrix (genes x samples) -> Normalization -> DESeq2/edgeR -> Volcano/MA -> GSEA/ORA Unit Formula Use case RPKM/FPKM (C / L) / N x 1e9 Single/paired-end, within-sample only TPM (C/L) / sum(Cj/Lj) x 1e6 Cross-sample comparison (sums to 1M) DESeq2 size factor median(count / gene_geomean) per sample Robust normalization for DE testing C = read count, L = gene length (bp), N = total mapped reads. Never feed TPM/RPKM/FPKM into DESeq2/edgeR — use raw counts; report TPM only for visualization. Goal: get a per-sample scale factor that corrects for library size and composition without being skewed by a few very highly expressed genes. Approach: DESeq2's median-of-ratios — divide each gene's count by its across-sample geometric mean, then take the per-sample median of those ratios. import numpy as np import pandas as pd def deseq2_size_factors ( count_matrix: pd.DataFrame ) -> pd.Series: """Compute DESeq2-style median-of-ratios size factors. count_matrix: genes (rows) x samples (columns), raw integer counts. Only genes with nonzero counts in every sample contribute to the geometric mean (required since log(0) is undefined). """ nonzero_mask = (count_matrix > 0 ). all (axis= 1 ) filtered = count_matrix.loc[nonzero_mask] geo_means = np.exp(np.log(filtered).mean(axis= 1 )) return pd.Series( {s: np.median(filtered[s] / geo_means) for s in filtered.columns} ) def counts_to_tpm ( counts: np.ndarray, lengths: np.ndarray ) -> np.ndarray: """Convert raw counts to TPM (length-normalize, then scale to 1e6).""" rate = counts / lengths return rate / rate. sum () * 1e6 size_factors = deseq2_size_factors(counts_df) normalized_counts = counts_df.div(size_factors, axis= 1 ) Goal: call differentially expressed genes with correct FDR control. Approach: for production use DESeq2 (R) or pydeseq2 (Python) — both fit a negative-binomial GLM per gene and borrow strength across genes for dispersion. A quick non-parametric fallback (Mann-Whitney + BH) is fine for prototyping only. library ( DESeq2 ) dds <- DESeqDataSetFromMatrix ( countData = count_matrix , colData = sample_info , design = ~ condition ) dds <- DESeq ( dds ) # size factors + dispersion + Wald test res <- results ( dds , contrast = c ( "condition" , "Treatment" , "Control" ) ) res_shrunk <- lfcShrink ( dds , coef = "condition_Treatment_vs_Control" , type = "apeglm" ) sig <- subset ( res_shrunk , padj < 0.05 & abs ( log2FoldChange ) > 1 ) summary ( res ) from pydeseq2.dds import DeseqDataSet from pydeseq2.ds import DeseqStats dds = DeseqDataSet(counts=count_matrix, metadata=sample_info, design_factors= "condition" ) dds.deseq2() stat_res = DeseqStats(dds, contrast=[ "condition" , "Treatment" , "Control" ]) stat_res.summary() results_df = stat_res.results_df # baseMean, log2FoldChange, pvalue, padj from scipy import stats from statsmodels.stats.multitest import multipletests def simple_de ( count_df: pd.DataFrame, ctrl_samples, treat_samples ) -> pd.DataFrame: """Prototype-only DE test: BH-corrected Mann-Whitney on DESeq2-normalized counts. Underpowered vs. DESeq2/edgeR (no dispersion borrowing) — do not use for publication. """ sf = deseq2_size_factors(count_df) norm = count_df.div(sf, axis= 1 ) rows = [] for gene in count_df.index: c, t = norm.loc[gene, ctrl_samples], norm.loc[gene, treat_samples] lfc = np.log2((t.mean() + 1 ) / (c.mean() + 1 )) _, p = stats.mannwhitneyu(c, t, alternative= "two-sided" ) rows.append({ "gene" : gene, "log2FC" : lfc, "pvalue" : p, "baseMean" : (c.mean() + t.mean()) / 2 }) df = pd.DataFrame(rows).set_index( "gene" ) df[ "padj" ] = multipletests(df[ "pvalue" ], method= "fdr_bh" )[ 1 ] return df.sort_values( "pvalue" ) Goal: get from raw reads to a count matrix. Approach: align + count (splice-aware, gene-level counts) or pseudoalign + import (fast, transcript-level then aggregated). STAR + featureCounts (alignment-based): 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 *.bam Salmon (alignment-free; aggregate transcript counts to genes with tximport in R): 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 Pitfalls TPM/RPKM/FPKM for DE testing — ratio-of-ratios artifacts and no variance model; always pass raw counts to DESeq2/edgeR. Library composition bias — one highly-expressed gene can suppress the apparent expression of everything else; DESeq2 median-of-ratios and edgeR TMM are robust, plain CPM/TPM are not. Skipping dispersion shrinkage — DESeq2/edgeR borrow information across genes to stabilize low-count variance estimates; a per-gene t-test or Mann-Whitney misses this and is underpowered, especially at n<3/group. Too few replicates — use >=3 biological replicates per condition (5+ preferred); replicates improve power far more than deeper sequencing. No LFC shrinkage — raw log2FoldChange is noisy for low-count genes; use lfcShrink(..., type="apeglm") before ranking/plotting. Not filtering low counts — genes with near-zero counts across samples inflate the multiple-testing burden and destabilize dispersion estimates; filter before testing. Confusing biological with technical replicates — technical replicates (same sample, resequenced) understate true variance and inflate significance. See Also bio-differential-expression-deseq2-basics — DESeq2 API details and design formulas bio-differential-expression-edger-basics — edgeR/TMM alternative workflow bio-rna-quantification-tximport-workflow — Salmon/kallisto to gene-level counts bio-pathway-analysis-gsea — downstream gene set enrichment on DE results