Transcriptomics

Definition (RNA, transcriptome). RNA (ribonucleic acid) is a single-stranded polymer over the alphabet $\Sigma_{RNA} = \{\texttt{A}, \texttt{C}, \texttt{G}, \texttt{U}\}$, where uracil (U) replaces thymine (T). RNA is synthesized by RNA polymerase using DNA as a template in a process called transcription. The RNA sequence is identical to the coding (sense) strand of the DNA (with T replaced by U). The transcriptome of a cell at time $t$ is the multiset of all RNA molecules present in the cell at $t$.

Definition (RNA classes). The major classes of RNA in eukaryotes are:

• mRNA (messenger RNA): encodes a protein. It carries information from DNA to the ribosome. In eukaryotes it has a 5' cap, a coding sequence (CDS), and a 3' poly-A tail.
• rRNA (ribosomal RNA): structural and catalytic component of ribosomes. The most abundant RNA by mass ($\sim$80% of total RNA).
• tRNA (transfer RNA): adaptor molecule that delivers amino acids during translation. Each tRNA has an anticodon complementary to one codon.
• ncRNA (non-coding RNA): a broad class including lncRNA (long non-coding RNA, $>$200 nt), miRNA (microRNA, $\sim$22 nt, post-transcriptional repressors), and snRNA (small nuclear RNA, involved in splicing).

Definition (Pre-mRNA, introns, exons, splicing). In eukaryotes, a gene is not a contiguous coding sequence. The primary transcript (pre-mRNA) contains alternating exons (retained in the mature mRNA) and introns (removed). The process of removing introns and joining exons is called splicing, catalyzed by the spliceosome. Let a gene have exons $E_1, \ldots, E_k$ of lengths $\ell_1, \ldots, \ell_k$. The mature mRNA has length $\sum_{i=1}^k \ell_i$, which can be much shorter than the genomic span of the gene (sometimes by a factor of 10–100 in humans, where introns average $\sim$3,400 bp vs. exons of $\sim$170 bp).

complement = str.maketrans('ACGT', 'TGCA')

def transcribe(dna_coding_strand):
    return dna_coding_strand.replace('T', 'U')

def splice(pre_mrna, intron_intervals):
    # intron_intervals: list of (start, end) 0-based half-open
    exons, pos = [], 0
    for start, end in sorted(intron_intervals):
        exons.append(pre_mrna[pos:start])
        pos = end
    exons.append(pre_mrna[pos:])
    return ''.join(exons)

Definition (Paired-end reads, insert size). In paired-end sequencing, both ends of each cDNA fragment are sequenced, producing two reads (called mates) per fragment. If the forward read starts at position $i$ and the reverse read ends at position $j$ on the transcript, the insert size is $j - i + 1$. Paired-end information (1) confirms that both mates align to the same gene or isoform, (2) spans exon-exon junctions when individual reads are too short, and (3) enables estimation of the insert size distribution. In single-end sequencing, only one read is produced per fragment, at lower cost but with less information.

Definition (Stranded library). In a stranded (strand-specific) library, the protocol preserves which DNA strand the RNA was transcribed from. The first-strand cDNA is synthesized with dUTP incorporated into the second strand, which is then degraded before PCR, so only the first strand (same orientation as the original RNA) is amplified. Each read's strand of origin is therefore known. This disambiguates overlapping genes on opposite strands and improves quantification of antisense transcription. In an unstranded library, strand information is lost.

Definition (Sequencing depth for RNA-seq). The sequencing depth $N$ is the total number of reads in a sample. Unlike genome sequencing, optimal depth in RNA-seq depends on the biological question. For standard differential expression of abundant genes: $N \approx 20\text{--}30 \times 10^6$ reads. For detection of lowly expressed genes or rare isoforms: $N \geq 100 \times 10^6$. The critical constraint is not coverage uniformity but count precision: a gene expressed at 0.01% of total mRNA receives only $\sim$100 reads at $N = 10^6$, giving high Poisson noise.

from collections import defaultdict

def build_kmer_index(transcripts, k):
    # transcripts: dict {transcript_id: sequence}
    index = defaultdict(set)
    for tid, seq in transcripts.items():
        for i in range(len(seq) - k + 1):
            index[seq[i:i+k]].add(tid)
    return index

def pseudoalign(read, index, k):
    compatible = None
    for i in range(len(read) - k + 1):
        hits = index.get(read[i:i+k])
        if hits:
            compatible = hits if compatible is None else compatible & hits
    return compatible or set()

Definition (Raw counts, count matrix). After alignment, reads are assigned to genes (or transcripts). The raw count $c_{gj}$ is the number of reads in sample $j$ assigned to gene $g$. The count matrix $C \in \mathbb{N}^{G \times J}$ collects all raw counts across $G$ genes and $J$ samples. Raw counts cannot be directly compared between genes (longer genes attract more reads) or between samples (different library sizes produce different total counts). Normalization is required before any comparison.

Definition (RPKM, FPKM, TPM). Let $c_g$ be the raw count of gene $g$, $\ell_g$ its effective length in kilobases, and $N$ the total number of mapped reads in millions. Three normalization units are in common use:

• RPKM (Reads Per Kilobase per Million): $\text{RPKM}_g = c_g / (\ell_g \cdot N)$. Corrects for gene length and library size.
• FPKM (Fragments Per Kilobase per Million): same as RPKM but counts fragments (pairs) rather than individual reads. Used for paired-end data.
• TPM (Transcripts Per Million): $\text{TPM}_g = \frac{c_g / \ell_g}{\sum_{g'} c_{g'} / \ell_{g'}} \times 10^6$. TPM values sum to exactly $10^6$ per sample, making them directly comparable between samples. TPM is now preferred over RPKM/FPKM.

def compute_tpm(counts, lengths_kb):
    rates = [c / l for c, l in zip(counts, lengths_kb)]
    total = sum(rates)
    return [r / total * 1e6 for r in rates]

def compute_rpkm(counts, lengths_kb, total_mapped_millions):
    return [c / (l * total_mapped_millions) for c, l in zip(counts, lengths_kb)]

Definition (Negative binomial model). The negative binomial distribution $\text{NB}(\mu, \alpha)$ with mean $\mu > 0$ and dispersion $\alpha \geq 0$ has variance $\text{Var}(X) = \mu + \alpha \mu^2$. When $\alpha = 0$ it reduces to Poisson. The extra term $\alpha \mu^2$ captures overdispersion: at higher mean expression, absolute variance grows quadratically. In the DESeq2 model, the count $c_{gj}$ for gene $g$ in sample $j$ satisfies $$c_{gj} \sim \text{NB}(\mu_{gj},\, \alpha_g), \qquad \mu_{gj} = s_j \cdot q_{gj}$$ where $s_j$ is a sample-specific size factor (accounting for different library sizes) and $q_{gj}$ is the normalized expression level. The gene-specific dispersion $\alpha_g$ is shared across all samples of the same gene.

Definition (Size factors, median-of-ratios). The size factor $s_j$ of sample $j$ corrects for sequencing depth and RNA composition differences. DESeq2 estimates it as follows. For each gene $g$, compute the geometric mean across samples: $\bar{c}_g = \bigl(\prod_{j=1}^J c_{gj}\bigr)^{1/J}$. The size factor is the median of per-gene ratios: $$s_j = \operatorname{median}_g \left\{\frac{c_{gj}}{\bar{c}_g}\right\}$$ Genes with $\bar{c}_g = 0$ are excluded. Dividing raw counts by $s_j$ gives the normalized counts $\tilde{c}_{gj} = c_{gj} / s_j$, comparable across samples. This estimator is robust to highly expressed outlier genes (which would bias simpler depth normalization).

import math
from statistics import median

def size_factors(count_matrix):
    # count_matrix: list of rows, each row = counts for one gene across all samples
    n_samples = len(count_matrix[0])
    sf = []
    for j in range(n_samples):
        ratios = []
        for row in count_matrix:
            # geometric mean over samples for this gene
            nonzero = [c for c in row if c > 0]
            if len(nonzero) == n_samples:
                gm = math.exp(sum(math.log(c) for c in row) / n_samples)
                ratios.append(row[j] / gm)
        sf.append(median(ratios))
    return sf

Definition (Log fold-change, Wald test). For gene $g$, the log$_2$ fold-change (LFC) between condition $B$ and condition $A$ is $\text{LFC}_g = \log_2(\hat{q}_{gB} / \hat{q}_{gA})$, where $\hat{q}_{gA}$ and $\hat{q}_{gB}$ are estimated mean normalized expressions. A positive LFC means higher expression in $B$. The Wald test tests $H_0: \text{LFC}_g = 0$ by forming $$W_g = \frac{\widehat{\text{LFC}}_g}{\widehat{\text{SE}}(\widehat{\text{LFC}}_g)}$$ which under $H_0$ is approximately standard Normal. The two-sided $p$-value is $p_g = 2\Phi(-|W_g|)$, where $\Phi$ is the standard Normal CDF.

Definition (Multiple testing, BH correction). Testing $G \approx 20{,}000$ genes simultaneously inflates false positives: even with no true signal, $\sim$1,000 genes would have $p < 0.05$ by chance. The false discovery rate (FDR) is $\text{FDR} = \mathbb{E}[V/R]$, where $V$ is the number of false rejections and $R$ the total number of rejections. The Benjamini--Hochberg (BH) procedure controls FDR at level $q$: sort $p$-values $p_{(1)} \leq \cdots \leq p_{(G)}$; let $k^* = \max\{k : p_{(k)} \leq kq/G\}$; reject all hypotheses $1, \ldots, k^*$. The adjusted $p$-value of gene $(k)$ is $\tilde{p}_{(k)} = \min_{j \geq k}(G \cdot p_{(j)} / j)$.

def bh_correction(pvalues, q=0.05):
    n = len(pvalues)
    indexed = sorted(enumerate(pvalues), key=lambda x: x[1])
    adjusted = [0.0] * n
    min_so_far = 1.0
    for rank, (i, p) in reversed(list(enumerate(indexed, 1))):
        adj = min(p * n / rank, 1.0)
        min_so_far = min(min_so_far, adj)
        adjusted[i] = min_so_far
    return adjusted

Definition (Splice graph). Given splice-aligned reads on a genomic locus, the splice graph $SG = (V, E)$ is a directed acyclic graph where each node $v \in V$ represents a maximal contiguous genomic interval covered by reads without interruption (an exonic segment), and each directed edge $(u, v) \in E$ represents an exon-exon junction observed in at least one spliced read. Each path from a source node to a sink node in $SG$ corresponds to a candidate transcript isoform. The assembly problem reduces to: which paths through $SG$ best explain the observed read coverage?

def em_isoform(compatible_sets, n_isoforms, n_iter=100):
    # compatible_sets[r] = set of isoform indices compatible with read r
    theta = [1.0 / n_isoforms] * n_isoforms
    for _ in range(n_iter):
        counts = [0.0] * n_isoforms
        for S in compatible_sets:
            denom = sum(theta[k] for k in S)
            for k in S:
                counts[k] += theta[k] / denom
        total = sum(counts)
        theta = [c / total for c in counts]
    return theta

Definition (UMI count matrix). After demultiplexing by cell barcode and UMI deduplication, the output is a UMI count matrix $X \in \mathbb{N}^{G \times C}$ where $X_{gc}$ is the number of distinct UMIs observed for gene $g$ in cell $c$. The matrix is extremely sparse: a typical cell detects $\sim$2,000--5,000 genes out of $\sim$30,000, so $\sim$85--93\% of entries are zero. This sparsity has two causes: genuine absence of expression, and technical dropout (a transcript was present but not captured).

Definition (Log-normalization). Each cell's raw counts are divided by its library size $N_c = \sum_g X_{gc}$ and scaled by a factor $S = 10{,}000$ (CP10K -- counts per 10,000): $\tilde{X}_{gc} = (X_{gc} / N_c) \times S$. A log-transformation is applied: $Y_{gc} = \log_2(\tilde{X}_{gc} + 1)$. The $+1$ pseudocount prevents $\log(0)$ and stabilizes variance for low counts. The resulting matrix $Y \in \mathbb{R}^{G \times C}$ is the starting point for all downstream analyses.

import math
from statistics import mean, variance

def log_normalize(count_matrix, scale=1e4):
    # count_matrix[c][g] = UMI count for cell c, gene g
    result = []
    for row in count_matrix:
        total = sum(row)
        result.append([math.log2(x / total * scale + 1) for x in row])
    return result

def highly_variable_genes(log_norm, top_n=2000):
    n_cells = len(log_norm)
    n_genes = len(log_norm[0])
    col = [[log_norm[c][g] for c in range(n_cells)] for g in range(n_genes)]
    var_per_gene = [variance(col[g]) if n_cells > 1 else 0 for g in range(n_genes)]
    ranked = sorted(range(n_genes), key=lambda g: -var_per_gene[g])
    return set(ranked[:top_n])

Definition (kNN graph and clustering). After PCA, each cell is a point in $\mathbb{R}^k$. A $k$-nearest-neighbor graph (kNN graph) is constructed: each cell is connected to its $k$ nearest neighbors in PCA space (typically $k = 15$, Euclidean distance). Community detection on this graph identifies clusters of cells with similar transcriptional profiles. The Leiden algorithm optimizes a modularity-like quality function over graph partitions. A resolution parameter $\gamma > 0$ controls granularity: larger $\gamma$ yields more, smaller clusters.

Definition (Gene Ontology). The Gene Ontology (GO) is a controlled vocabulary describing gene functions, structured as a directed acyclic graph. Each GO term is a node representing a biological concept; edges represent "is-a" or "part-of" relationships. GO has three independent ontologies: Biological Process (BP, e.g., "DNA repair"), Molecular Function (MF, e.g., "ATP binding"), and Cellular Component (CC, e.g., "nucleus"). The true path rule: if gene $g$ is annotated with term $t$, it is implicitly annotated with all ancestors of $t$ in the DAG.

Definition (Over-representation analysis, hypergeometric test). Given a hit list $S$ (e.g., differentially expressed genes) of size $n$ drawn from a universe of $N$ genes, and a GO term $T$ annotating $K$ genes in the universe, the overlap $k = |S \cap T|$ follows a hypergeometric distribution under the null that $S$ is a random subset: $$P(X = k) = \frac{\binom{K}{k}\binom{N-K}{n-k}}{\binom{N}{n}}$$ The one-sided $p$-value for over-representation is $P(X \geq k) = \sum_{x=k}^{\min(K,n)} P(X = x)$. After BH correction over all tested terms, enriched terms with small adjusted $p$-value represent biological processes disproportionately represented in the hit list.

from math import comb

def hypergeometric_pvalue(N, K, n, k):
    # P(X >= k), X ~ Hypergeometric(N, K, n)
    return sum(
        comb(K, x) * comb(N - K, n - x) / comb(N, n)
        for x in range(k, min(K, n) + 1)
    )

def over_representation(hit_set, gene_sets, universe):
    N = len(universe)
    n = len(hit_set & universe)
    results = []
    for term, genes in gene_sets.items():
        K = len(genes & universe)
        k = len(hit_set & genes & universe)
        if K > 0 and k > 0:
            p = hypergeometric_pvalue(N, K, n, k)
            results.append((term, k, K, p))
    return sorted(results, key=lambda x: x[3])

Definition (Gene Set Enrichment Analysis, GSEA). ORA requires a hard threshold to define the hit list. GSEA avoids this: it takes a ranked list of all $G$ genes (ranked by signed LFC or $-\log_{10}(p)$) and a gene set $T$. Walking down the ranked list, an enrichment score ES is computed as a running sum that increments by $1/|T|$ when a gene in $T$ is encountered and decrements by $1/(G - |T|)$ otherwise. ES($T$) is the maximum deviation from zero. A large positive ES means genes in $T$ cluster at the top of the list (preferentially up-regulated). Significance is assessed by permutation of gene labels.

def gsea_score(ranked_genes, gene_set):
    n = len(ranked_genes)
    n_set = sum(1 for g in ranked_genes if g in gene_set)
    n_miss = n - n_set
    if n_set == 0 or n_miss == 0:
        return 0.0
    running, max_dev = 0.0, 0.0
    for g in ranked_genes:
        running += 1.0 / n_set if g in gene_set else -1.0 / n_miss
        if abs(running) > abs(max_dev):
            max_dev = running
    return max_dev