# Bio Applied Rna Seq Analysis

> Bulk RNA-seq — STAR/HISAT2/featureCounts or Salmon/kallisto quantification, TPM/DESeq2 size-factor normalization, DESeq2/pydeseq2 DE testing. Use for RNA-seq design, count matrices, or DE analysis with DESeq2, edgeR, pydeseq2.

- Skill: `pavel-kravchenko/bio-applied-rna-seq-analysis` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-applied-rna-seq-analysis`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-applied-rna-seq-analysis/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Coding & Dev Tools
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/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

```text
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 |

```bash
# 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.

```python
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

```python
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.

```r
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)
```

```python
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