# Core bioinformatics & visualisation libraries
 # install.packages(c("ggplot2","dplyr","tidyr","scales","knitr","kableExtra","gridExtra","RColorBrewer","pheatmap","ggrepel"))
 # BiocManager::install(c("DESeq2","edgeR","apeglm","limma","GenomicRanges","Biostrings","VariantAnnotation","clusterProfiler","org.Hs.eg.db"))
suppressPackageStartupMessages({
  library(ggplot2)
  library(dplyr)
  library(tidyr)
  library(scales)
  library(knitr)
  library(kableExtra)
  library(RColorBrewer)
  library(gridExtra)
  library(ggrepel)
  library(apeglm)
})

How to use this document:
Every step includes: (1) Historical context — when it was first achieved and by whom; (2) Plain-English explanation + underlying algorithm/statistics; (3) Exact file formats with annotated examples; (4) R and Python code a bioinformatician would actually run; (5) Industry interview questions from top genomics companies; and (6) Full citations with hyperlinks.
Command syntax follows each tool’s official manual — always check --help or current docs before production use, as flags evolve.

1 PART 0 — Historical Milestones: How We Got Here

1.1 The Road to Modern Genomics

Understanding why we do things the way we do requires knowing how we got here. The sequencing revolution did not happen overnight.

1.2 Timeline of Key Historical Achievements

milestones <- data.frame(
  Year = c(1953, 1977, 1977, 1986, 1990, 1995, 2000, 2001, 2003, 2005, 2007, 2008,
           2009, 2010, 2012, 2013, 2014, 2016, 2018, 2019, 2020, 2022, 2023, 2024),
  Event = c(
    "Watson & Crick: DNA double helix (Nature)",
    "Sanger: dideoxy chain-termination sequencing (PNAS)",
    "Maxam & Gilbert: chemical cleavage sequencing",
    "Leroy Hood: automated fluorescent Sanger sequencer (ABI 370A)",
    "Human Genome Project launched (NIH/DOE)",
    "First complete bacterial genome: H. influenzae (Science)",
    "Celera Genomics draft human genome (Craig Venter)",
    "HGP + Celera papers in Nature & Science; reference genome published",
    "Human Genome Project completed; 99.99% accuracy",
    "454 pyrosequencing: first 'next-generation' platform commercially available",
    "Illumina GA: short-read SBS; 1000 Genomes Project launched",
    "BWA aligner published; SOLiD sequencing",
    "SAMtools + SAM/BAM format; GATK v1; RNA-seq method established",
    "FastQC released; 1000 Genomes pilot paper",
    "STAR aligner; Seurat v1; Oxford Nanopore Technology founded",
    "ENCODE Project data released; DESeq2 published",
    "Illumina HiSeq X Ten: $1000 genome achieved; ClinVar launched",
    "DeepMind AlphaFold v1; 10x Genomics Chromium scRNA-seq",
    "DeepVariant; GATK4; gnomAD v2",
    "PacBio HiFi (CCS) published; Cellranger 3.0",
    "AlphaFold2 (Nature); COVID-19 rapid sequencing achieves SARS-CoV-2 genome in 10h",
    "T2T-CHM13: first truly complete human genome (Science); Nanopore R10.4",
    "AlphaFold database covers >200M proteins; gnomAD v4 (800K genomes)",
    "Illumina NovaSeq X Plus; long-read clinical sequencing FDA cleared"
  ),
  Category = c("Biology","Sequencing","Sequencing","Technology","Initiative","Sequencing",
                "Initiative","Initiative","Initiative","Technology","Technology","Algorithms",
                "Algorithms","Tools","Tools","Analysis","Technology","Technology","Tools",
                "Technology","AI/Tools","Genomics","Databases","Technology"),
  Importance = c(10,10,8,9,10,9,9,10,10,9,9,8,9,7,9,8,10,9,8,9,10,10,9,9)
)

cat_colors <- c(
  "Biology"="#e74c3c","Sequencing"="#3498db","Technology"="#2ecc71",
  "Initiative"="#f39c12","Algorithms"="#9b59b6","Tools"="#1abc9c",
  "Analysis"="#e67e22","Databases"="#34495e","AI/Tools"="#e91e63"
)

ggplot(milestones, aes(x = Year, y = Importance, color = Category, label = substr(Event,1,45))) +
  geom_segment(aes(x=Year, xend=Year, y=0, yend=Importance), linewidth=0.5, alpha=0.4) +
  geom_point(aes(size=Importance), alpha=0.85) +
  geom_text_repel(size=2.5, max.overlaps=20, segment.size=0.3, box.padding=0.3) +
  scale_color_manual(values=cat_colors) +
  scale_size_continuous(range=c(3,8)) +
  scale_x_continuous(breaks=seq(1950,2025,5)) +
  labs(title="Key Milestones in Sequencing & Computational Genomics (1953–2024)",
       x="Year", y="Historical Importance (subjective scale)",
       color="Category", caption="Sources: primary papers cited throughout this document") +
  theme_minimal(base_size=11) +
  theme(legend.position="bottom", axis.text.x=element_text(angle=45,hjust=1),
        plot.title=element_text(face="bold"))
Major milestones in sequencing and computational genomics history

Major milestones in sequencing and computational genomics history

1.2.1 Key Historical Highlights

1977 — Sanger Sequencing Revolution
Frederick Sanger (Cambridge) and Walter Gilbert (Harvard) independently invented sequencing methods in 1977. Sanger’s dideoxy (chain-termination) method dominated for 30 years and sequenced the human genome. It worked by using fluorescently labelled dideoxy nucleotides (ddNTPs) that lacked the 3′-OH needed to extend the chain, terminating synthesis at each base; separation by gel electrophoresis then read the ladder of fragments by size.
Reference: Sanger F., Nicklen S., Coulson A.R. (1977). “DNA sequencing with chain-terminating inhibitors.” PNAS 74(12):5463–5467. [Link]

1986 — Automated Sequencing
Leroy Hood and colleagues at Caltech coupled fluorescent dyes to Sanger chemistry and used a laser detector, enabling the first automated sequencer (ABI 370A). This reduced sequencing from a manual, week-long process to a semi-automated overnight run.
Reference: Smith L.M. et al. (1986). “Fluorescence detection in automated DNA sequence analysis.” Nature 321:674–679.

1990–2003 — Human Genome Project (HGP)
The $3 billion, 13-year international effort sequenced all 3.2 billion base pairs of the human genome using Sanger sequencing on an industrial scale (~1,500 Sanger sequencers running simultaneously). The final reference (GRCh37/hg19) was declared complete in April 2003. Key achievement: establishing the reference that all modern variant calling compares against.
Reference: Lander E.S. et al. (2001). “Initial sequencing and analysis of the human genome.” Nature 409:860–921.

2005–2007 — The NGS Revolution
454 Life Sciences (Roche) launched the first commercial next-generation sequencer in 2004/2005 using pyrosequencing (detection of pyrophosphate released when a nucleotide is incorporated). Illumina’s acquisition of Solexa in 2007 and the launch of the Genome Analyzer brought bridge-amplification + sequencing-by-synthesis to market, eventually reducing the cost of sequencing a human genome from $10M (2007) to $1,000 (2014) to <$200 (2024).
Reference: Margulies M. et al. (2005). “Genome sequencing in microfabricated high-density picolitre reactors.” Nature 437:376–380.

2022 — The Complete Human Genome (T2T)
The “complete” HGP reference had ~8% of the genome missing — primarily telomeric, centromeric, and highly repetitive regions that short reads couldn’t assemble. The Telomere-to-Telomere (T2T) Consortium used PacBio HiFi and Oxford Nanopore long reads to finally assemble these regions, producing the truly complete CHM13 reference in 2022.
Reference: Nurk S. et al. (2022). “The complete sequence of a human genome.” Science 376:44–53.

2 PART 1 — Sequencing Platforms: Mechanisms, Outputs, and Use Cases

2.1 Sequencing Platforms Overview

platforms <- data.frame(
  Platform = c("Illumina NovaSeq X Plus", "Oxford Nanopore PromethION", "PacBio Revio (HiFi)",
                "10x Genomics Chromium", "MGI DNBSEQ-T7", "Element Bio AVITI",
                "Ultima Genomics UG 100"),
  Mechanism = c(
    "Bridge amplification on flow cell; SBS with reversible terminators; 4-colour fluorescence per cycle",
    "Single molecule through protein nanopore; ionic current disruption decoded by neural-net basecaller",
    "Circular Consensus Sequencing (CCS): SMRT cell ZMW wells; same circularised fragment read 10–30×",
    "GEM droplets: 10x barcode bead + single cell; pooled Illumina sequencing of barcoded cDNA",
    "DNA Nanoball (DNB) rolling circle amplification on array; combinatorial Probe-Anchor Binding (cPAL)",
    "Avidity sequencing: avidite reagents bind multiple polymerases on same cluster; improved accuracy",
    "Flow cell with periodic dNTPs; rolling circle amplification; $1 per genome claimed"
  ),
  Read_Length = c("100–300 bp, PE", "N50 > 50 kb routinely, up to Mb", "15–25 kb, ~Q30",
                   "28 bp barcode + 90 bp cDNA R2", "50–300 bp, PE", "150 bp, PE", "300–500 bp"),
  Raw_Output = c("~10 Tb/run (flow cell)", "~290 Gb/flow cell (24 cells/run)", "~90 Gb/SMRT cell (8/run)",
                  "~100M read pairs/lane", "~7 Tb/run", "~800 Gb/run", "~1 Tb/run"),
  Use_Case = c("WGS, WES, RNA-seq, scRNA-seq, amplicons", "SV, de novo assembly, direct RNA, metagenomics",
               "Long-read WGS, full-length isoforms (Iso-Seq), methylation", "scRNA-seq, scATAC-seq, spatial",
               "WGS, WES, NIPT", "WGS, WES, RNA-seq", "Population WGS at ultra-low cost"),
  stringsAsFactors = FALSE
)
kable(platforms, caption="Major sequencing platforms (2024)", booktabs=TRUE) %>%
  kable_styling(bootstrap_options=c("striped","hover","condensed"), full_width=TRUE, font_size=12) %>%
  column_spec(1, bold=TRUE, color="white", background="#2980b9")
Major sequencing platforms (2024)
Platform Mechanism Read_Length Raw_Output Use_Case
Illumina NovaSeq X Plus Bridge amplification on flow cell; SBS with reversible terminators; 4-colour fluorescence per cycle 100–300 bp, PE ~10 Tb/run (flow cell) WGS, WES, RNA-seq, scRNA-seq, amplicons
Oxford Nanopore PromethION Single molecule through protein nanopore; ionic current disruption decoded by neural-net basecaller N50 > 50 kb routinely, up to Mb ~290 Gb/flow cell (24 cells/run) SV, de novo assembly, direct RNA, metagenomics
PacBio Revio (HiFi) Circular Consensus Sequencing (CCS): SMRT cell ZMW wells; same circularised fragment read 10–30× 15–25 kb, ~Q30 ~90 Gb/SMRT cell (8/run) Long-read WGS, full-length isoforms (Iso-Seq), methylation
10x Genomics Chromium GEM droplets: 10x barcode bead + single cell; pooled Illumina sequencing of barcoded cDNA 28 bp barcode + 90 bp cDNA R2 ~100M read pairs/lane scRNA-seq, scATAC-seq, spatial
MGI DNBSEQ-T7 DNA Nanoball (DNB) rolling circle amplification on array; combinatorial Probe-Anchor Binding (cPAL) 50–300 bp, PE ~7 Tb/run WGS, WES, NIPT
Element Bio AVITI Avidity sequencing: avidite reagents bind multiple polymerases on same cluster; improved accuracy 150 bp, PE ~800 Gb/run WGS, WES, RNA-seq
Ultima Genomics UG 100 Flow cell with periodic dNTPs; rolling circle amplification; $1 per genome claimed 300–500 bp ~1 Tb/run Population WGS at ultra-low cost

2.2 Phred Quality Score — Deep Dive

The Phred quality score is the foundational unit of trust in sequencing data.

\[Q = -10 \times \log_{10}(P_{error})\]

Q Score P(error) Accuracy Interpretation
Q10 1 in 10 90% Unacceptable for most analyses
Q20 1 in 100 99% Minimum acceptable threshold
Q30 1 in 1,000 99.9% Standard Illumina benchmark
Q40 1 in 10,000 99.99% Excellent; PacBio HiFi target
Q50 1 in 100,000 99.999% Near-perfect; clinical standards

Origin: Ewing B. & Green P. (1998). “Base-calling of automated sequencer traces using phred II.” Genome Research 8:186–194. [PubMed]

# Phred score encoding demonstration
phred_to_prob <- function(q) 10^(-q/10)
q_scores <- seq(0, 41, 1)
probs <- phred_to_prob(q_scores)
ascii_chars <- sapply(q_scores, function(q) rawToChar(as.raw(q + 33)))

phred_df <- data.frame(
  Q = q_scores,
  P_error = probs,
  ASCII_char = ascii_chars,
  Accuracy_pct = (1 - probs) * 100
)

ggplot(phred_df, aes(x=Q, y=Accuracy_pct)) +
  geom_line(color="#2980b9", linewidth=1.5) +
  geom_point(data=phred_df %>% filter(Q %in% c(10,20,30,40)), 
             aes(color=factor(Q)), size=4) +
  geom_text_repel(data=phred_df %>% filter(Q %in% c(10,20,30,40)),
                  aes(label=paste0("Q",Q,"\n",round(Accuracy_pct,2),"%")),
                  size=3.5, fontface="bold") +
  scale_y_continuous(limits=c(85,100.1), labels=function(x) paste0(x,"%")) +
  scale_color_brewer(palette="Set1", guide="none") +
  labs(title="Phred Quality Score vs. Base Accuracy",
       subtitle="Q30 (1-in-1000 error rate) is the standard Illumina benchmark",
       x="Phred Quality Score (Q)", y="Base Call Accuracy (%)") +
  theme_minimal(base_size=12)

# ASCII encoding table
kable(head(phred_df, 15), digits=6, caption="Phred+33 ASCII encoding (first 15 Q values)") %>%
  kable_styling(bootstrap_options="striped", full_width=FALSE)
Phred+33 ASCII encoding (first 15 Q values)
Q P_error ASCII_char Accuracy_pct
0 1.000000 ! 0.00000
1 0.794328 ” 20.56718
2 0.630957 # 36.90427
3 0.501187 $ 49.88128
4 0.398107 % 60.18928
5 0.316228 & 68.37722
6 0.251189 ’ 74.88114
7 0.199526 ( 80.04738
8 0.158489 ) 84.15107
9 0.125893
87.41075
10 0.100000
90.00000
11 0.079433 , 92.05672
12 0.063096
93.69043
13 0.050119 . 94.98813
14 0.039811 / 96.01893 
# Python equivalent — Phred decoding from a FASTQ quality string
def decode_phred(qual_string, offset=33):
    """Decode a FASTQ quality string into Phred scores and error probabilities."""
    q_scores = [ord(c) - offset for c in qual_string]
    p_errors = [10**(-q/10) for q in q_scores]
    return list(zip(qual_string, q_scores, [round(p,6) for p in p_errors]))

qual = "#FFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFF:"
decoded = decode_phred(qual)
print(f"{'Char':^5} {'Q':^5} {'P(error)':^10}")
print("-" * 25)
for char, q, p in decoded[:8]:
    print(f"{char:^5} {q:^5} {p:^10.6f}")
# '#' = Q2 (very low, often position 1 of Illumina read — known artifact)
# 'F' = Q37, ':'= Q25

3 PART 1A — FASTQ Format: Complete Anatomy

3.1 FASTQ Structure

@A00123:45:HG7NKDSXX:1:1101:5678:1000 1:N:0:CGATGT
NAGCTGACGTTTGCAAGGCTAGCATGCATGCATGCATG
+
#FFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFF:
Line Prefix Content Detail
1 @ Read identifier Instrument:Run:FlowCell:Lane:Tile:X:Y SPACE Mate:Filter:Control:Index
2 — Base sequence ACGTN; N = undetermined base
3 + Separator Optionally repeats header (legacy practice)
4 — Quality string ASCII-encoded Phred+33; one character per base of line 2

Paired-end files: sample_R1.fastq.gz and sample_R2.fastq.gz — read n in R1 and read n in R2 are the two ends of the same DNA fragment. The insert size (fragment length between read starts) is critical to alignment confidence.
Reference: Cock P.J.A. et al. (2010). “The Sanger FASTQ file format.” Nucleic Acids Research 38(6):1767–1771. [Link]

# R: parse and visualise quality scores from a FASTQ quality string
qual_string <- "#FFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFF:"
q_scores <- utf8ToInt(qual_string) - 33
p_errors <- 10^(-q_scores/10)

qc_df <- data.frame(
  Position = seq_along(q_scores),
  Q = q_scores,
  P_error = p_errors,
  Pass = q_scores >= 30
)

ggplot(qc_df, aes(x=Position, y=Q, fill=Pass)) +
  geom_bar(stat="identity", width=0.8) +
  geom_hline(yintercept=30, linetype="dashed", color="red", linewidth=1) +
  annotate("text", x=5, y=31.5, label="Q30 threshold", color="red", size=3.5) +
  scale_fill_manual(values=c("TRUE"="#27ae60","FALSE"="#e74c3c"),
                    labels=c("TRUE"="Pass Q30","FALSE"="Fail Q30")) +
  labs(title="Per-base Quality Scores for Example Read",
       x="Position in Read", y="Phred Quality Score", fill="") +
  theme_minimal(base_size=12)

# Python: full FASTQ parser with statistics
import gzip
from collections import defaultdict
import statistics

def parse_fastq(filepath, max_reads=1000):
    """Parse FASTQ file, return per-position quality stats."""
    opener = gzip.open if filepath.endswith('.gz') else open
    pos_quals = defaultdict(list)
    n_reads = 0
    with opener(filepath, 'rt') as fh:
        while True:
            header = fh.readline().strip()
            if not header: break
            seq = fh.readline().strip()
            fh.readline()  # '+'
            qual = fh.readline().strip()
            for i, c in enumerate(qual):
                pos_quals[i].append(ord(c) - 33)
            n_reads += 1
            if n_reads >= max_reads: break
    
    stats = {}
    for pos, qs in sorted(pos_quals.items()):
        stats[pos] = {
            'mean': statistics.mean(qs),
            'median': statistics.median(qs),
            'q10': sorted(qs)[int(0.10*len(qs))],
            'q90': sorted(qs)[int(0.90*len(qs))]
        }
    return stats, n_reads

# Usage:
# stats, n = parse_fastq("sample_R1.fastq.gz")
# print(f"Parsed {n} reads; median Q at pos 0: {stats[0]['median']}")

3.2 Quality Control (QC)

3.2.1 FastQC — Metrics Every Bioinformatician Must Know

fastqc sample_R1.fastq.gz sample_R2.fastq.gz \
  --threads 8 \
  --outdir qc/ \
  --extract

Industry Q: “Walk me through a FastQC report — what do you look for first?” (Illumina, Novartis, Genentech)

  1. Per-base sequence quality — should be ≥Q30 across all positions; quality drop at 3’ end is normal
  2. Per-tile sequence quality (Illumina) — hot spots indicate flow cell lane issues
  3. Per-sequence GC content — bimodal distribution suggests contamination or adapter dimers
  4. Sequence duplication levels — >50% indicates low library complexity (wet-lab problem)
  5. Adapter content — if >10% at any position, trimming is mandatory
  6. Overrepresented sequences — BLAST these; often adapter sequences or rRNA contamination

Reference: Andrews S. (2010). FastQC: A Quality Control tool. Babraham Bioinformatics. [Link]

# Simulate FastQC-style per-base quality plot
set.seed(42)
n_positions <- 150
pos <- 1:n_positions
mean_q <- c(rep(38, 10), seq(38, 36, length.out=40), rep(36, 60), seq(36, 28, length.out=40))
q10  <- mean_q - 8 + rnorm(n_positions, 0, 1)
q25  <- mean_q - 4 + rnorm(n_positions, 0, 0.5)
q75  <- mean_q + 2 + rnorm(n_positions, 0, 0.5)
q90  <- mean_q + 5 + rnorm(n_positions, 0, 1)

fqc_df <- data.frame(pos=pos, mean=mean_q, q10=q10, q25=q25, q75=q75, q90=q90)

ggplot(fqc_df, aes(x=pos)) +
  annotate("rect", xmin=0, xmax=n_positions+1, ymin=28, ymax=30, fill="#ffd700", alpha=0.3) +
  annotate("rect", xmin=0, xmax=n_positions+1, ymin=20, ymax=28, fill="#ff6b6b", alpha=0.2) +
  annotate("rect", xmin=0, xmax=n_positions+1, ymin=30, ymax=42, fill="#90ee90", alpha=0.2) +
  geom_ribbon(aes(ymin=q10, ymax=q90), fill="#3498db", alpha=0.2) +
  geom_ribbon(aes(ymin=q25, ymax=q75), fill="#3498db", alpha=0.4) +
  geom_line(aes(y=mean), color="#e74c3c", linewidth=1.2) +
  geom_hline(yintercept=30, linetype="dashed", color="grey30") +
  annotate("text", x=10, y=40.5, label="Very Good (Q>30)", size=3, color="darkgreen") +
  annotate("text", x=10, y=29, label="Acceptable (Q28-30)", size=3, color="goldenrod4") +
  annotate("text", x=10, y=24, label="Poor (Q<28)", size=3, color="red") +
  scale_y_continuous(limits=c(18,43), breaks=seq(20,40,5)) +
  labs(title="FastQC-Style Per-Base Quality Score Distribution (Simulated 150bp PE Read)",
       subtitle="Red line = mean quality; blue band = 10th–90th percentile",
       x="Position in Read (bp)", y="Phred Quality Score") +
  theme_minimal(base_size=12)

3.2.2 MultiQC — Aggregating Across Samples

multiqc . \
  --outdir multiqc_report/ \
  --title "Project XYZ QC Summary" \
  --filename multiqc_report.html \
  --ignore-samples "negative_control*"

Reference: Ewels P. et al. (2016). “MultiQC: summarize analysis results for multiple tools and samples in a single report.” Bioinformatics 32(19):3047–3048. [Link]

3.2.3 Trimming — fastp

fastp \
  -i sample_R1.fastq.gz -I sample_R2.fastq.gz \
  -o trim_R1.fastq.gz   -O trim_R2.fastq.gz \
  --detect_adapter_for_pe \
  --qualified_quality_phred 20 \
  --unqualified_percent_limit 40 \
  --length_required 36 \
  --cut_right --cut_window_size 4 --cut_mean_quality 20 \
  --thread 16 \
  -j fastp.json -h fastp.html

Interview Q (Illumina/10x Genomics): “What is the difference between quality trimming and adapter trimming, and when is each necessary?”

Answer: Adapter trimming removes synthetic adapter sequences that were ligated to library ends but read through when the insert is shorter than the read length (critical for small RNA-seq, ATAC-seq). Quality trimming removes low-quality bases from 3′ ends (common in all data). For WGS with inserts >150bp, adapter contamination is minimal but quality trimming still improves mapping. For small RNA-seq with 20–22nt inserts, adapter trimming is essential. fastp’s --detect_adapter_for_pe uses overlap analysis to auto-detect adapters — the overlap between R1 and R2 reveals how much is adapter vs. insert.

Reference: Chen S. et al. (2018). “fastp: an ultra-fast all-in-one FASTQ preprocessor.” Bioinformatics 34(17):i884–i890. [Link]

Alternative trimming tools:

Tool Approach Best For
Trim Galore Cutadapt wrapper Bisulfite (WGBS)
Trimmomatic Sliding window Legacy pipelines
fastp Overlap detection General-purpose (fastest)
BBDuk k-mer-based Contamination screening
Cutadapt Regex adapter matching Custom adapters

4 PART 2 — DNA Alignment: From FASTQ to Analysis-Ready BAM

4.1 Reference Genome Preparation

refs <- data.frame(
  Build = c("GRCh38/hg38","GRCh37/hg19","T2T-CHM13v2.0","GRCh38 + ALT contigs","CHM13 + Y chromosome"),
  Released = c("2013","2009","2022","2019","2023"),
  Size_Gb = c("3.1","3.1","3.1","3.2","3.1"),
  Completeness = c("~92% (gaps at centromeres/telomeres)","~91%","100% (first truly complete)","92%+ALT haplotypes","100%"),
  Recommended_Use = c("Standard clinical/research WGS","Legacy; still common in clinical","Research; SV/repeat studies","Population diversity studies","Most complete reference"),
  stringsAsFactors = FALSE
)
kable(refs, caption="Human reference genome builds") %>%
  kable_styling(bootstrap_options=c("striped","hover"), full_width=TRUE)
Human reference genome builds
Build Released Size_Gb Completeness Recommended_Use
GRCh38/hg38 2013 3.1 ~92% (gaps at centromeres/telomeres) Standard clinical/research WGS
GRCh37/hg19 2009 3.1 ~91% Legacy; still common in clinical
T2T-CHM13v2.0 2022 3.1 100% (first truly complete) Research; SV/repeat studies
GRCh38 + ALT contigs 2019 3.2 92%+ALT haplotypes Population diversity studies
CHM13 + Y chromosome 2023 3.1 100% Most complete reference 
# --- Reference genome indexing (one-time setup) ---

# BWA-MEM2 index (faster than BWA-MEM for large genomes)
bwa-mem2 index GRCh38.fa
# Produces: .amb, .ann, .bwt.2bit.64, .pac, .0123

# SAMtools FASTA index
samtools faidx GRCh38.fa
# Produces: GRCh38.fa.fai — tab-delimited: name, length, byte_offset, bases_per_line, bytes_per_line

# GATK sequence dictionary
gatk CreateSequenceDictionary -R GRCh38.fa
# Produces: GRCh38.dict — SAM-format @SQ header lines for all chromosomes

# Verify index integrity
samtools view -H sorted.bam | grep "^@SQ" | wc -l  # should match reference chromosome count

Reference genome: Genome Reference Consortium GRCh38.p14 [NCBI]; Nurk S. et al. (2022). Science 376:44–53. [Link]

4.2 Alignment Algorithms — Deep Dive

4.2.1 Burrows-Wheeler Transform (BWT) — the Mathematics

Why BWA is fast: the FM-index (based on the Burrows-Wheeler Transform) allows searching a 3-billion-base genome for a 150-bp read in O(read_length) time — essentially independent of genome size.

The BWT algorithm (1994): 1. Create all cyclic rotations of the genome string T (length n)
2. Sort them lexicographically → produces the Suffix Array (SA)
3. The last column of the sorted matrix (BWT) has a special property: runs of identical bases cluster together → enables RLE compression (same principle as bzip2)
4. The rank (occ) and select data structures on BWT enable backward search: narrow the SA interval one character at a time, finding all exact matches in O(m) where m = pattern length

Reference: Burrows M., Wheeler D.J. (1994). “A block-sorting lossless data compression algorithm.” DEC Technical Report 124. [PDF]

# Python: BWT and backward search demo (educational — shows the O(m) pattern)
def bwt(s):
    """Compute Burrows-Wheeler Transform of string s (append '$' as sentinel)."""
    s = s + '$'
    rotations = sorted([s[i:] + s[:i] for i in range(len(s))])
    return ''.join(r[-1] for r in rotations), rotations

def backward_search(bwt_str, pattern, first_col_chars):
    """BWT backward search: find all occurrences of pattern in O(m)."""
    # Simplified demo using linear rank — real FM-index uses wavelet trees
    top, bot = 0, len(bwt_str)
    for c in reversed(pattern):
        top = sum(1 for x in bwt_str[:top] if x == c) + first_col_chars.get(c, 0)
        bot = sum(1 for x in bwt_str[:bot] if x == c) + first_col_chars.get(c, 0)
        if top >= bot:
            return 0  # pattern not found
    return bot - top  # number of occurrences

genome = "ACGTACGTACGT"
bwt_str, rotations = bwt(genome)
print(f"Genome:  {genome}")
print(f"BWT:     {bwt_str}")

# Count first-column characters (needed for backward search)
sorted_bwt = sorted(bwt_str)
first_col = {}
for i, c in enumerate(sorted_bwt):
    if c not in first_col:
        first_col[c] = i

hits = backward_search(bwt_str, "ACG", first_col)
print(f"Occurrences of 'ACG': {hits}")  # Should find 3 in "ACGTACGTACGT"

4.2.2 BWA-MEM2 Alignment

# --- Alignment with BWA-MEM2 (2–3× faster than BWA-MEM, identical results) ---
bwa-mem2 mem \
  -t 16 \
  -R "@RG\tID:sample1_L001\tSM:sample1\tLB:library1\tPL:ILLUMINA\tPU:A00123.45.HG7NKDSXX.1" \
  -K 100000000 \
  -Y \
  GRCh38.fa \
  trim_R1.fastq.gz trim_R2.fastq.gz \
  | samtools sort -@ 8 -m 2G -o sample1.sorted.bam -

samtools index -@ 8 sample1.sorted.bam

# Validate the BAM header
samtools view -H sample1.sorted.bam | head -20

Critical: the Read Group (-R) flag explained

ID = unique lane identifier; SM = sample name (used by GATK for joint calling); LB = library (MarkDuplicates operates per-library to correctly handle re-sequenced libraries); PL = platform (BQSR uses this); PU = platform unit (flowcell.lane for per-lane error correction in BQSR).

Missing -R breaks: MarkDuplicates (can’t deduplicate per-library), BQSR (per-sample error models), HaplotypeCaller multi-sample calling. Fixing after the fact requires re-alignment — extremely expensive at WGS scale.
Reference: GATK Read Groups documentation. [Link]

For long reads (Nanopore/PacBio):

# PacBio HiFi
minimap2 -ax map-hifi \
  -t 16 \
  -R "@RG\tID:pb_sample\tSM:sample1\tPL:PACBIO" \
  GRCh38.fa \
  reads.hifi.fastq.gz \
  | samtools sort -@ 8 -o pb_sorted.bam -

# Oxford Nanopore (R10.4, Q20 chemistry)
minimap2 -ax map-ont \
  -t 16 \
  --secondary=no \
  GRCh38.fa \
  nanopore_reads.fastq.gz \
  | samtools sort -@ 8 -o ont_sorted.bam -

References: Li H. (2013). arXiv:1303.3997 [BWA-MEM]; Vasimuddin Md. et al. (2019). “Efficient Architecture-Aware Acceleration of BWA-MEM for Multicore Systems.” IPDPS 2019 [BWA-MEM2]; Li H. (2018). Bioinformatics 34(18):3094–3100. [minimap2]

4.3 SAM/BAM Format — Complete Anatomy

flags <- data.frame(
  Bit = c(1,2,4,8,16,32,64,128,256,512,1024,2048),
  Hex = c("0x1","0x2","0x4","0x8","0x10","0x20","0x40","0x80","0x100","0x200","0x400","0x800"),
  Description = c(
    "Read is paired (PE library)",
    "Read is in a proper pair (both mapped, correct orientation & distance)",
    "Read is unmapped",
    "Mate is unmapped",
    "Read is on reverse strand",
    "Mate is on reverse strand",
    "Read is read 1 (first in pair)",
    "Read is read 2 (second in pair)",
    "Alignment is a secondary alignment (alternative alignment position)",
    "Read failed quality controls (not passing filter)",
    "Read is a PCR or optical duplicate (MarkDuplicates sets this)",
    "Alignment is a supplementary alignment (chimeric read)"
  ),
  Common_Values = c("Always set for PE","Expect >95% for WGS","Check for unmapped reads","","Filter with -F 16","","Set for R1","Set for R2","Exclude with -F 256","Exclude with -F 512","Exclude with -F 1024","Usually exclude with -F 2048")
)
kable(flags, caption="SAM FLAG field bitmask — every bioinformatician must know these", booktabs=TRUE) %>%
  kable_styling(bootstrap_options=c("striped","hover","condensed"), full_width=TRUE, font_size=12) %>%
  column_spec(1, bold=TRUE)
SAM FLAG field bitmask — every bioinformatician must know these
Bit Hex Description Common_Values
1 0x1 Read is paired (PE library) Always set for PE
2 0x2 Read is in a proper pair (both mapped, correct orientation & distance) Expect >95% for WGS
4 0x4 Read is unmapped Check for unmapped reads
8 0x8 Mate is unmapped
16 0x10 Read is on reverse strand Filter with -F 16
32 0x20 Mate is on reverse strand
64 0x40 Read is read 1 (first in pair) Set for R1
128 0x80 Read is read 2 (second in pair) Set for R2
256 0x100 Alignment is a secondary alignment (alternative alignment position) Exclude with -F 256
512 0x200 Read failed quality controls (not passing filter) Exclude with -F 512
1024 0x400 Read is a PCR or optical duplicate (MarkDuplicates sets this) Exclude with -F 1024
2048 0x800 Alignment is a supplementary alignment (chimeric read) Usually exclude with -F 2048 
# R function to decode SAM flags
decode_sam_flag <- function(flag) {
  flag_meanings <- c(
    "1"="read paired", "2"="proper pair", "4"="read unmapped",
    "8"="mate unmapped", "16"="read reverse strand", "32"="mate reverse strand",
    "64"="read1", "128"="read2", "256"="secondary", "512"="qcfail",
    "1024"="duplicate", "2048"="supplementary"
  )
  set_bits <- names(which(sapply(as.integer(names(flag_meanings)),
                                  function(b) bitwAnd(flag, b) > 0)))
  cat(sprintf("FLAG %d decodes to: %s\n", flag,
              paste(flag_meanings[set_bits], collapse=", ")))
}

decode_sam_flag(99)   # 1+2+32+64 — paired, proper pair, mate reverse, read1
#> FLAG 99 decodes to:
decode_sam_flag(147)  # 1+2+16+128 — paired, proper pair, read reverse, read2
#> FLAG 147 decodes to:
decode_sam_flag(1024) # duplicate only
#> FLAG 1024 decodes to:
# Python: parse a BAM file and compute coverage statistics using pysam
import pysam
import collections

def bam_qc_metrics(bam_path, region=None):
    """
    Compute key QC metrics from a BAM file.
    Equivalent to: samtools flagstat + samtools idxstats
    """
    bam = pysam.AlignmentFile(bam_path, "rb")
    
    metrics = collections.Counter()
    insert_sizes = []
    
    for read in bam.fetch(region=region):
        metrics['total'] += 1
        if read.is_paired:            metrics['paired'] += 1
        if read.is_proper_pair:       metrics['proper_pair'] += 1
        if read.is_unmapped:          metrics['unmapped'] += 1
        if read.is_duplicate:         metrics['duplicate'] += 1
        if read.is_secondary:         metrics['secondary'] += 1
        if read.is_supplementary:     metrics['supplementary'] += 1
        if read.mapping_quality >= 20: metrics['mapq_ge20'] += 1
        
        # Collect insert sizes (only for read 1, proper pairs)
        if (read.is_read1 and read.is_proper_pair and 
            not read.is_unmapped and not read.mate_is_unmapped):
            if 0 < abs(read.template_length) < 1000:
                insert_sizes.append(abs(read.template_length))
    
    bam.close()
    
    print("\n=== BAM QC Metrics ===")
    for k, v in metrics.items():
        pct = 100*v/metrics['total'] if metrics['total'] > 0 else 0
        print(f"  {k:20s}: {v:>10,d} ({pct:.1f}%)")
    
    if insert_sizes:
        import statistics
        print(f"\n  Insert size (n={len(insert_sizes):,}):")
        print(f"    Median: {statistics.median(insert_sizes):.0f} bp")
        print(f"    Mean:   {statistics.mean(insert_sizes):.0f} bp")
    
    return metrics, insert_sizes

# Usage:
# metrics, inserts = bam_qc_metrics("analysis_ready.bam")

4.4 Mark Duplicates

Why duplicates arise: PCR amplifies the same original DNA fragment before sequencing. Both copies map identically (same chromosome, start position, end position, strand). They are not independent evidence for a variant — counting them would artificially inflate confidence.

MarkDuplicates logic:
1. Groups reads sharing identical 5′ genomic coordinates (per-library, using the LB read-group tag)
2. Within each duplicate set, keeps the read with the highest base-quality sum as the “representative”
3. Sets FLAG bit 1024 on all others — does not delete them (still auditable)
4. Optical duplicates (clusters too physically close on the flow cell, from a single original cluster being miscalled as two) are also flagged separately

# MarkDuplicatesSpark (GATK4 — parallelized, preferred for large files)
gatk MarkDuplicatesSpark \
  -I sorted.bam \
  -O dedup.bam \
  -M dedup_metrics.txt \
  --spark-master local[16] \
  --optical-duplicate-pixel-distance 2500   # 2500 for patterned flow cells (NovaSeq), 100 for non-patterned

# Check metrics
cat dedup_metrics.txt | grep -A 2 "ESTIMATED_LIBRARY_SIZE"
# Typical duplication rates by library type and input amount
dup_data <- data.frame(
  Library_Type = rep(c("WGS 30x (1ug input)", "WGS 30x (50ng input)", "WES 100x",
                        "RNA-seq (low input)", "scRNA-seq 10x", "ATAC-seq"), 1),
  Typical_Dup_Rate = c(5, 35, 20, 40, 60, 50),
  Warning_Threshold = c(20, 50, 35, 60, 75, 65),
  stringsAsFactors = FALSE
)

ggplot(dup_data, aes(x=reorder(Library_Type, Typical_Dup_Rate), y=Typical_Dup_Rate)) +
  geom_bar(stat="identity", aes(fill=Typical_Dup_Rate > 25), width=0.6) +
  geom_point(aes(y=Warning_Threshold), shape=23, size=4, fill="orange", color="darkred") +
  coord_flip() +
  scale_fill_manual(values=c("FALSE"="#27ae60","TRUE"="#e74c3c"),
                    labels=c("FALSE"="Acceptable","TRUE"="High — check input"), guide=FALSE) +
  labs(title="Typical PCR Duplication Rates by Library Type",
       subtitle="Orange diamond = warning threshold; bars = typical observed rate",
       x="", y="Duplication Rate (%)") +
  theme_minimal(base_size=12)

Reference: Picard toolkit, Broad Institute. [Link]
Van der Auwera G.A. & O’Connor B.D. (2020). Genomics in the Cloud. O’Reilly.

4.5 Base Quality Score Recalibration (BQSR) — Why and How

The problem BQSR solves:
Illumina sequencers have systematic, reproducible error biases — e.g., the dinucleotide “GpC” context is systematically harder to call correctly; quality scores at cycle 74 of a 75-cycle run are slightly worse than reported. These biases are consistent across reads but not captured in the raw Phred score. If not corrected, they inflate false-positive variant calls.

BQSR algorithm:
1. BaseRecalibrator: compares observed mismatches against a “known variation” VCF (dbSNP + Mills indels) — mismatches at known-variant sites are excluded from the error model (they’re real variants, not errors). Builds a recalibration table stratified by read group, reported quality score, sequence context (2 bp before/after), and machine cycle.
2. ApplyBQSR: rewrites quality scores in the BAM using the correction table.
3. AnalyzeCovariates (optional): generates before/after calibration plots for QC.

# Step 1: Build recalibration model
gatk BaseRecalibrator \
  -I dedup.bam \
  -R GRCh38.fa \
  --known-sites dbsnp_146.hg38.vcf.gz \
  --known-sites Mills_and_1000G_gold_standard.indels.hg38.vcf.gz \
  --known-sites 1000G_phase1.snps.high_confidence.hg38.vcf.gz \
  -O recal.table

# Step 2: Apply recalibration
gatk ApplyBQSR \
  -I dedup.bam \
  -R GRCh38.fa \
  --bqsr-recal-file recal.table \
  -O analysis_ready.bam

# Step 3 (Optional): Visualize recalibration
gatk AnalyzeCovariates \
  -before recal.table \
  -after recal_after.table \
  -plots BQSR_plots.pdf
# Visualize BQSR effect: reported Q vs empirical Q before and after
set.seed(123)
reported_q <- 10:40
empirical_before <- reported_q + rnorm(31, -3, 1.5)  # systematic underestimation
empirical_after  <- reported_q + rnorm(31, 0, 0.4)   # near-perfect calibration

bqsr_df <- data.frame(
  Reported = rep(reported_q, 2),
  Empirical = c(empirical_before, empirical_after),
  Stage = rep(c("Before BQSR","After BQSR"), each=31)
)

ggplot(bqsr_df, aes(x=Reported, y=Empirical, color=Stage)) +
  geom_abline(slope=1, intercept=0, linetype="dashed", color="grey50", linewidth=1) +
  geom_line(linewidth=1.5) +
  geom_point(size=2) +
  scale_color_manual(values=c("Before BQSR"="#e74c3c","After BQSR"="#27ae60")) +
  labs(title="BQSR Effect: Reported vs. Empirical Quality Scores",
       subtitle="Dashed line = perfect calibration; Before BQSR shows systematic underestimation",
       x="Reported Phred Quality Score", y="Empirical Quality Score", color="") +
  theme_minimal(base_size=12) +
  theme(legend.position="top")

Reference: DePristo M.A. et al. (2011). “A framework for variation discovery and genotyping using next-generation DNA sequencing data.” Nature Genetics 43:491–498. [Link]

4.6 BAM QC — What Numbers Matter

bam_qc <- data.frame(
  Metric = c("Mean coverage depth (WGS)", "Mean coverage depth (WES)",
             "% properly paired reads", "% mapped reads", "Duplication rate",
             "Insert size median (PE)", "% on-target (WES)", "% bases ≥ Q30",
             "Ti/Tv ratio (WGS)", "Ti/Tv ratio (WES)"),
  Acceptable = c("≥25×","≥50×","≥95%","≥95%","<30%","150–500bp","≥70%","≥75%","2.0–2.1","3.0–3.3"),
  Concerning = c("<20×","<30×","<90%","<85%",">40%","<100 or >700bp","<60%","<70%","<1.8 or >2.3","<2.5 or >3.8"),
  Tool = c("mosdepth/samtools","mosdepth","samtools flagstat","samtools flagstat",
           "Picard MarkDuplicates","Picard CollectInsertSizeMetrics",
           "Picard HsMetrics","FastQC","GATK/vcftools","GATK/vcftools"),
  stringsAsFactors=FALSE
)
kable(bam_qc, caption="BAM/variant QC thresholds for WGS and WES") %>%
  kable_styling(bootstrap_options=c("striped","hover"), full_width=TRUE) %>%
  column_spec(3, color="white", background="#e74c3c") %>%
  column_spec(2, color="white", background="#27ae60")
BAM/variant QC thresholds for WGS and WES
Metric Acceptable Concerning Tool
Mean coverage depth (WGS) ≥25× <20× mosdepth/samtools
Mean coverage depth (WES) ≥50× <30× mosdepth
% properly paired reads ≥95% <90% samtools flagstat
% mapped reads ≥95% <85% samtools flagstat
Duplication rate <30% >40% Picard MarkDuplicates
Insert size median (PE) 150–500bp <100 or >700bp Picard CollectInsertSizeMetrics
% on-target (WES) ≥70% <60% Picard HsMetrics
% bases ≥ Q30 ≥75% <70% FastQC
Ti/Tv ratio (WGS) 2.0–2.1 <1.8 or >2.3 GATK/vcftools
Ti/Tv ratio (WES) 3.0–3.3 <2.5 or >3.8 GATK/vcftools 
# Comprehensive BAM QC commands
# 1. Coverage depth (fastest tool)
mosdepth --threads 8 --quantize 0:1:10:50: sample_cov analysis_ready.bam

# 2. Flagstat
samtools flagstat -@ 8 analysis_ready.bam > flagstat.txt

# 3. Insert size distribution
picard CollectInsertSizeMetrics \
  INPUT=analysis_ready.bam \
  OUTPUT=insert_size_metrics.txt \
  HISTOGRAM_FILE=insert_size_histogram.pdf

# 4. WES target coverage
picard CollectHsMetrics \
  INPUT=analysis_ready.bam \
  BAIT_INTERVALS=targets.interval_list \
  TARGET_INTERVALS=targets.interval_list \
  OUTPUT=hs_metrics.txt \
  REFERENCE_SEQUENCE=GRCh38.fa

# 5. Per-sample summary (fast)
samtools stats analysis_ready.bam | grep "^SN" | cut -f 2,3 | head -30

5 PART 3 — Variant Calling: SNVs, Indels, SVs, and CNVs

5.1 Variant Calling Overview

var_types <- data.frame(
  Variant_Class = c("SNV (Single Nucleotide Variant)", "Small Indel (<50bp)", 
                     "Large Deletion (>50bp)", "Duplication/Amplification",
                     "Inversion", "Translocation/BND", "Copy Number Variant (CNV)",
                     "Mobile Element Insertion", "Repeat Expansion"),
  Size = c("1 bp","1–49 bp",">50 bp",">50 bp",">50 bp","Chromosomal","Kb–Mb","100–7000 bp","Variable"),
  Typical_Frequency = c("~4M/genome","~600K/genome","~2,500/genome","~1,000/genome",
                         "~1,000/genome","~700/genome","~500/genome","~1,000/genome","~30 pathogenic/genome"),
  Tools = c("GATK HaplotypeCaller, DeepVariant, Strelka2",
             "GATK HaplotypeCaller, DeepVariant",
             "Manta, PBSV, Sniffles2 (long read)",
             "Manta, CNVkit, PURPLE",
             "Manta, PBSV, SVABA",
             "DELLY, LUMPY, PBSV",
             "CNVkit, GATK gCNV, Canvas",
             "MELT, TLDR",
             "ExpansionHunter, Straglr, TRGT"),
  stringsAsFactors=FALSE
)
kable(var_types, caption="Variant types, sizes, and detection tools") %>%
  kable_styling(bootstrap_options=c("striped","hover","condensed"), full_width=TRUE, font_size=12)
Variant types, sizes, and detection tools
Variant_Class Size Typical_Frequency Tools
SNV (Single Nucleotide Variant) 1 bp ~4M/genome GATK HaplotypeCaller, DeepVariant, Strelka2
Small Indel (<50bp) 1–49 bp ~600K/genome GATK HaplotypeCaller, DeepVariant
Large Deletion (>50bp) >50 bp ~2,500/genome Manta, PBSV, Sniffles2 (long read)
Duplication/Amplification >50 bp ~1,000/genome Manta, CNVkit, PURPLE
Inversion >50 bp ~1,000/genome Manta, PBSV, SVABA
Translocation/BND Chromosomal ~700/genome DELLY, LUMPY, PBSV
Copy Number Variant (CNV) Kb–Mb ~500/genome CNVkit, GATK gCNV, Canvas
Mobile Element Insertion 100–7000 bp ~1,000/genome MELT, TLDR
Repeat Expansion Variable ~30 pathogenic/genome ExpansionHunter, Straglr, TRGT

5.2 GATK HaplotypeCaller — Algorithm Deep Dive

HaplotypeCaller’s key innovation over pileup callers:

Rather than examining one position at a time, HaplotypeCaller:
1. Identifies “active regions” — windows with >20% of reads showing evidence of variation
2. Builds a De Bruijn graph of all k-mer paths through those reads
3. Enumerates candidate haplotypes by walking graph paths
4. Re-aligns every read against every candidate haplotype using a Pair Hidden Markov Model (PairHMM)
5. Applies Bayes’ theorem to compute posterior genotype probabilities

This local assembly correctly resolves nearby indels that confuse pileup callers.

Reference: Poplin R. et al. (2018). “Scaling accurate genetic variant discovery to tens of thousands of samples.” bioRxiv 201178. [Link]

# ---- Germline variant calling ----

# Single sample → GVCF mode (required for joint genotyping)
gatk HaplotypeCaller \
  -R GRCh38.fa \
  -I analysis_ready.bam \
  -O sample.g.vcf.gz \
  -ERC GVCF \
  --dbsnp dbsnp_146.hg38.vcf.gz \
  -L intervals.bed \
  --sample-name sample1 \
  --native-pair-hmm-threads 4

# ---- Joint genotyping across cohort ----
# Step 1: Import GVCFs into GenomicsDB (one per genomic interval for parallelism)
gatk GenomicsDBImport \
  --sample-name-map sample_map.txt \
  --genomicsdb-workspace-path /path/to/db/chr1 \
  -L chr1 \
  --reader-threads 5 \
  --batch-size 50

# Step 2: Joint genotyping
gatk GenotypeGVCFs \
  -R GRCh38.fa \
  -V gendb:///path/to/db/chr1 \
  -O cohort_chr1.vcf.gz \
  --dbsnp dbsnp_146.hg38.vcf.gz

# Step 3: Merge chromosomes
gatk GatherVcfs \
  $(ls cohort_chr*.vcf.gz | sed 's/^/-I /') \
  -O cohort_all.vcf.gz

Interview Q (GATK/Broad Institute style): “Why use GVCF mode + joint genotyping instead of calling each sample separately?”

Answer: Per-sample calling loses information at sites where a variant was below threshold. GVCF mode retains likelihood evidence at every site for every sample — even confident non-variant positions. Joint calling with GenomicsDBImport/GenotypeGVCFs can then distinguish true homozygous-reference (0/0 with deep, clean coverage) from a genuine missing call (insufficient coverage). This “statistical borrowing” across samples also improves sensitivity for variants that are sub-threshold in one sample but clearly called in others, and improves accuracy of rare variant genotyping in the cohort.

5.3 DeepVariant — AI-based Variant Calling

Historical context: DeepVariant (2018, Google Brain/Verily) was the first tool to apply deep learning (specifically ResNet-based image classification) to germline variant calling. Rather than a statistical model of sequencing errors, it treats the variant-calling problem as image classification: the pileup at each candidate site is rendered as an RGB image (3 channels: base identity, base quality, strand; rows = reads, columns = positions), and an Inception-v3 CNN trained on ~10M labeled examples from Genome in a Bottle truth sets classifies each site as hom-ref / het / hom-alt.
Reference: Poplin R. et al. (2018). “A universal SNP and small-indel variant caller using deep neural networks.” Nature Biotechnology 36:983–987. [Link]

# DeepVariant — runs as Docker/Singularity for reproducibility
BIN_VERSION="1.6.1"
docker run \
  -v "$(pwd):/input" \
  -v "$(pwd)/output:/output" \
  google/deepvariant:${BIN_VERSION} \
  /opt/deepvariant/bin/run_deepvariant \
  --model_type=WGS \
  --ref=/input/GRCh38.fa \
  --reads=/input/analysis_ready.bam \
  --output_vcf=/output/deepvariant.vcf.gz \
  --output_gvcf=/output/deepvariant.g.vcf.gz \
  --num_shards=32 \
  --intermediate_results_dir=/output/intermediate_results
# Performance comparison of variant callers (based on published benchmarks)
callers <- data.frame(
  Caller = rep(c("GATK HaplotypeCaller","DeepVariant","Strelka2","Clair3 (long-read)","PEPPER-Margin-DeepVariant"), 2),
  Variant_Type = rep(c("SNV","Indel"), each=5),
  F1_Score = c(0.9991, 0.9994, 0.9989, 0.9981, 0.9986,  # SNV F1
               0.9932, 0.9961, 0.9912, 0.9801, 0.9872),  # Indel F1
  stringsAsFactors = FALSE
)

ggplot(callers, aes(x=reorder(Caller, F1_Score), y=F1_Score, fill=Variant_Type)) +
  geom_bar(stat="identity", position="dodge", width=0.7) +
  geom_text(aes(label=round(F1_Score,4)), position=position_dodge(0.7), hjust=-0.1, size=3) +
  coord_flip() +
  scale_fill_manual(values=c("SNV"="#3498db","Indel"="#e74c3c")) +
  scale_y_continuous(limits=c(0.96,1.002), labels=function(x) sprintf("%.3f",x)) +
  labs(title="Variant Caller Performance (F1 Score) on GIAB HG002",
       subtitle="Based on published benchmarks; actual performance varies by data type and coverage",
       x="", y="F1 Score", fill="Variant Type",
       caption="Sources: GIAB benchmarks, individual tool papers") +
  theme_minimal(base_size=12) +
  theme(legend.position="top")

5.4 Somatic Variant Calling — Mutect2

History of somatic calling:
- 2013: MuTect (Cibulskis et al.) — first widely adopted somatic SNV caller using a Bayesian classifier
- 2019: Mutect2 (GATK4) — extends to indels, adds orientation bias filtering (critical for FFPE DNA), Panel of Normals approach
- Key challenge: somatic variants can be present in only a fraction of tumor cells (tumor heterogeneity) — may have variant allele fraction (VAF) as low as 1–5%

# ---- Somatic calling: tumor-normal pair ----

# Step 1: Create Panel of Normals (PoN) from unrelated healthy samples
gatk CreateSomaticPanelOfNormals \
  -V normal1.g.vcf.gz -V normal2.g.vcf.gz \
  -O pon.vcf.gz

# Step 2: Call somatic variants
gatk Mutect2 \
  -R GRCh38.fa \
  -I tumor.bam  -tumor tumor_sample_name \
  -I normal.bam -normal normal_sample_name \
  --panel-of-normals pon.vcf.gz \
  --germline-resource gnomad.hg38.vcf.gz \
  -O somatic_raw.vcf.gz \
  -bamout reassembled.bam

# Step 3: Filter artifacts
gatk GetPileupSummaries -I tumor.bam -V small_exac_common_3.hg38.vcf.gz -O tumor_pileups.table
gatk CalculateContamination -I tumor_pileups.table -O contamination.table

gatk FilterMutectCalls \
  -V somatic_raw.vcf.gz \
  --contamination-table contamination.table \
  --stats somatic_raw.vcf.gz.stats \
  -O somatic_filtered.vcf.gz

# Step 4: Orientation bias (critical for FFPE)
gatk LearnReadOrientationModel -I somatic_raw.vcf.gz.f1r2.tar.gz -O artifact_priors.tar.gz
# Then re-run FilterMutectCalls with --ob-priors artifact_priors.tar.gz

Reference: Cibulskis K. et al. (2013). Nature Biotechnology 31:213–219; Benjamin D. et al. (2019). bioRxiv 861054. [MuTect2]

5.5 Structural Variant Detection

# ---- Manta: short-read SV calling ----
configManta.py \
  --tumorBam tumor.bam \
  --normalBam normal.bam \
  --referenceFasta GRCh38.fa \
  --runDir manta_run/

python manta_run/runWorkflow.py -j 16 -g 32  # 16 threads, 32GB RAM

# Output: candidateSmallIndels.vcf.gz, candidateSV.vcf.gz, diploidSV.vcf.gz, somaticSV.vcf.gz

# ---- Sniffles2: long-read SV calling (Nanopore/PacBio) ----
sniffles \
  --input ont_sorted.bam \
  --vcf sniffles_sv.vcf.gz \
  --snf sniffles.snf \
  --threads 16 \
  --minsvlen 50 \
  --mapq 20

# ---- PBSV: PacBio-specific SV calling ----
pbsv discover --log-level INFO hifi.bam ref.svsig.gz
pbsv call GRCh38.fa ref.svsig.gz pbsv.vcf

# ---- CNVkit: copy number from targeted sequencing ----
cnvkit.py batch tumor.bam \
  --normal normal.bam \
  --targets targets.bed \
  --fasta GRCh38.fa \
  --access data/access-5k-mappable.hg38.bed \
  --output-reference reference.cnn \
  --output-dir cnvkit_output/ \
  -p 16

References: Chen X. et al. (2016). Bioinformatics 32(8):1220–1222; Talevich E. et al. (2016). PLOS Comp Biol 12(4):e1004873. [CNVkit]

5.6 VCF Format — Complete Anatomy

vcf_fields <- data.frame(
  Field = c("CHROM","POS","ID","REF","ALT","QUAL","FILTER","INFO","FORMAT","SAMPLE_GT"),
  Example = c("chr7","140453136","rs121913364","A","T","850.3","PASS",
               "DP=42;AF=0.50;AC=1;AN=2;DB;MQ=60",
               "GT:AD:DP:GQ:PL","0/1:20,22:42:99:880,0,710"),
  Description = c(
    "Chromosome (must match reference dict exactly)",
    "1-based position of REF allele start",
    "dbSNP/ClinVar rsID or '.' if novel",
    "Reference allele (must match reference genome at POS)",
    "Alternate allele(s) — comma-separated for multi-allelic sites",
    "Phred-scaled quality that the site is NOT a variant (higher = more confident variant)",
    "PASS or filter name(s) that failed; '.' if not filtered",
    "Semicolon-delimited key=value annotations: DP=depth, AF=allele freq, MQ=mapping quality",
    "Colon-delimited codes defining per-sample field order",
    "GT=0/1 het; AD=ref,alt depths; DP=total depth; GQ=genotype quality; PL=phred likelihoods"
  ),
  stringsAsFactors=FALSE
)
kable(vcf_fields, caption="VCF v4.2 field-by-field anatomy") %>%
  kable_styling(bootstrap_options=c("striped","hover","condensed"), full_width=TRUE, font_size=12) %>%
  column_spec(1, bold=TRUE, color="white", background="#2980b9") %>%
  column_spec(2, color="white", background="#34495e", monospace=TRUE)
VCF v4.2 field-by-field anatomy
Field Example Description
CHROM chr7 Chromosome (must match reference dict exactly)
POS 140453136 1-based position of REF allele start
ID rs121913364 dbSNP/ClinVar rsID or ‘.’ if novel
REF A Reference allele (must match reference genome at POS)
ALT T Alternate allele(s) — comma-separated for multi-allelic sites
QUAL 850.3 Phred-scaled quality that the site is NOT a variant (higher = more confident variant)
FILTER PASS PASS or filter name(s) that failed; ‘.’ if not filtered
INFO DP=42;AF=0.50;AC=1;AN=2;DB;MQ=60 Semicolon-delimited key=value annotations: DP=depth, AF=allele freq, MQ=mapping quality
FORMAT GT:AD:DP:GQ:PL Colon-delimited codes defining per-sample field order
SAMPLE_GT 0/1:20,22:42:99:880,0,710 GT=0/1 het; AD=ref,alt depths; DP=total depth; GQ=genotype quality; PL=phred likelihoods 
# Python: parse a VCF and compute Ti/Tv ratio using PyVCF or cyvcf2
import re

def compute_titv(vcf_path):
    """
    Compute Ti/Tv ratio from a VCF file.
    Transitions: A<->G (purines), C<->T (pyrimidines)
    Transversions: purine<->pyrimidine
    Expected: ~2.0-2.1 (WGS), ~3.0-3.3 (WES)
    """
    transitions = {frozenset(['A','G']), frozenset(['C','T'])}
    ti = tv = 0
    
    with open(vcf_path, 'r') as fh:
        for line in fh:
            if line.startswith('#'): continue
            fields = line.strip().split('\t')
            if len(fields) < 5: continue
            ref, alt = fields[3], fields[4]
            
            # Only SNVs (single nucleotide)
            if len(ref) != 1 or len(alt) != 1 or alt == '.': continue
            pair = frozenset([ref.upper(), alt.upper()])
            if pair in transitions: ti += 1
            elif len(pair) == 2:    tv += 1
    
    ratio = ti / tv if tv > 0 else float('inf')
    print(f"Transitions:    {ti:>8,}")
    print(f"Transversions:  {tv:>8,}")
    print(f"Ti/Tv ratio:    {ratio:.4f}")
    print(f"QC:             {'PASS' if 1.9 < ratio < 2.3 else 'WARN — check callset quality'}")
    return ratio

# Usage: ratio = compute_titv("filtered.vcf")
# Ti/Tv ratio across different filtering stringency levels (simulated)
titv_data <- data.frame(
  Filter_Level = c("No filter","QUAL>10","QUAL>30","VQSR 99%","VQSR 95%","Hard filter"),
  Ti_Tv = c(1.85, 1.95, 2.05, 2.09, 2.10, 2.08),
  N_Variants = c(5200000, 4800000, 4400000, 4200000, 4100000, 4350000)
)

p1 <- ggplot(titv_data, aes(x=reorder(Filter_Level, Ti_Tv), y=Ti_Tv)) +
  geom_bar(stat="identity", aes(fill=Ti_Tv > 2.0 & Ti_Tv < 2.2), width=0.6) +
  geom_hline(yintercept=c(2.0,2.1), linetype="dashed", color=c("orange","darkgreen")) +
  annotate("text", x=0.6, y=2.11, label="Ideal WGS", size=3, color="darkgreen") +
  coord_flip() +
  scale_fill_manual(values=c("FALSE"="#e74c3c","TRUE"="#27ae60"), guide=FALSE) +
  labs(title="Ti/Tv Ratio by Filter Stringency", x="", y="Ti/Tv Ratio") +
  theme_minimal(base_size=11)

p2 <- ggplot(titv_data, aes(x=reorder(Filter_Level, Ti_Tv), y=N_Variants/1e6)) +
  geom_bar(stat="identity", fill="#3498db", width=0.6) +
  coord_flip() +
  labs(title="Number of Variants Retained", x="", y="Variants (millions)") +
  theme_minimal(base_size=11)

grid.arrange(p1, p2, ncol=2)

5.7 Variant Filtering: Hard Filters vs VQSR vs CNN

# ---- Hard filtering (small cohorts, <30 samples) ----
# SNV filters
gatk VariantFiltration \
  -V cohort.vcf.gz \
  --filter-expression "QD < 2.0"          --filter-name "LowQD" \
  --filter-expression "FS > 60.0"         --filter-name "StrandBias" \
  --filter-expression "MQ < 40.0"         --filter-name "LowMQ" \
  --filter-expression "MQRankSum < -12.5" --filter-name "LowMQRS" \
  --filter-expression "ReadPosRankSum < -8.0" --filter-name "LowRPRS" \
  --filter-expression "SOR > 3.0"         --filter-name "HighSOR" \
  --select-type-to-include SNP \
  -O snv_filtered.vcf.gz

# ---- VQSR (large cohorts, >50 samples) ----
# Trains a Gaussian mixture model on known-truth sites
gatk VariantRecalibrator \
  -V cohort.vcf.gz \
  --resource:hapmap,known=false,training=true,truth=true,prior=15 hapmap.vcf.gz \
  --resource:omni,known=false,training=true,truth=false,prior=12 omni.vcf.gz \
  --resource:1000G,known=false,training=true,truth=false,prior=10 1000G_snps.vcf.gz \
  --resource:dbsnp,known=true,training=false,truth=false,prior=2 dbsnp.vcf.gz \
  -an QD -an MQRankSum -an ReadPosRankSum -an FS -an MQ -an SOR -an DP \
  -mode SNP \
  -O recal_snps.recal \
  --tranches-file recal_snps.tranches \
  --rscript-file recal_plots.R

gatk ApplyVQSR \
  -V cohort.vcf.gz \
  --recal-file recal_snps.recal \
  --tranches-file recal_snps.tranches \
  --truth-sensitivity-filter-level 99.5 \
  -mode SNP \
  -O vqsr_filtered.vcf.gz

6 PART 4 — Variant Annotation and Clinical Interpretation

6.1 VEP (Variant Effect Predictor) — Deep Dive

# Full VEP annotation with all major databases
vep \
  --input_file filtered.vcf.gz \
  --output_file annotated.vcf \
  --vcf \
  --cache --offline \
  --assembly GRCh38 \
  --everything \
  --fork 8 \
  --plugin CADD,GRCh38_CADD_v1.6.vcf.gz \
  --plugin SpliceAI,snv=spliceai_scores.masked.snv.hg38.vcf.gz \
  --plugin REVEL,revel_all_chromosomes.tsv.gz \
  --custom gnomad.genomes.v4.vcf.gz,gnomAD_v4,vcf,exact,0,AF,AF_afr,AF_eas,AF_nfe \
  --custom clinvar_20240301.vcf.gz,ClinVar,vcf,exact,0,CLNSIG,CLNDN \
  --fields "Uploaded_variation,Location,Allele,Gene,Feature,Consequence,IMPACT,gnomAD_AF,ClinVar_CLNSIG,CADD_PHRED,REVEL"

6.2 Pathogenicity Prediction Scores

scores_df <- data.frame(
  Score = c("CADD","REVEL","AlphaMissense","SIFT","PolyPhen-2","SpliceAI","SQUIRLS","EVE","PrimateAI"),
  Type = c("Combined","Missense","Missense","Missense","Missense","Splicing","Splicing","Missense","Missense"),
  Range = c("0–99 (Phred-scaled)","0–1","0–1 (AM score)","0–1 (lower=damaging)","0–1 (higher=damaging)",
             "0–1 per splice site","0–1","0–1","0–1"),
  Damaging_Threshold = c(">20 (top 1%),>30 (top 0.1%)","≥0.75","≥0.564 (likely pathogenic)","≤0.05","≥0.85",">0.5","≥0.5","≥0.5","≥0.803"),
  Method = c("Supervised ML + conservation","Ensemble of 18 tools","Protein LLM (AlphaFold-trained)","Sequence homology","Bayes naive classifier","Deep learning (splicing grammar)","Graph-based splicing","VAE evolutionary model","Deep CNN on primate variation"),
  Key_Reference = c("Kircher et al. 2014 NG","Ioannidis et al. 2016 Nature Genetics","Cheng et al. 2023 Science","Ng & Henikoff 2003 NAR","Adzhubei et al. 2010 NM","Jaganathan et al. 2019 Cell","Riepe et al. 2021 AJHG","Frazer et al. 2021 Nature","Sundaram et al. 2018 Nature Genetics"),
  stringsAsFactors=FALSE
)
kable(scores_df, caption="In silico pathogenicity prediction scores") %>%
  kable_styling(bootstrap_options=c("striped","hover","condensed"), full_width=TRUE, font_size=11) %>%
  column_spec(1, bold=TRUE)
In silico pathogenicity prediction scores
Score Type Range Damaging_Threshold Method Key_Reference
CADD Combined 0–99 (Phred-scaled) >20 (top 1%),>30 (top 0.1%) Supervised ML + conservation Kircher et al. 2014 NG
REVEL Missense 0–1 ≥0.75 Ensemble of 18 tools Ioannidis et al. 2016 Nature Genetics
AlphaMissense Missense 0–1 (AM score) ≥0.564 (likely pathogenic) Protein LLM (AlphaFold-trained) Cheng et al. 2023 Science
SIFT Missense 0–1 (lower=damaging) ≤0.05 Sequence homology Ng & Henikoff 2003 NAR
PolyPhen-2 Missense 0–1 (higher=damaging) ≥0.85 Bayes naive classifier Adzhubei et al. 2010 NM
SpliceAI Splicing 0–1 per splice site >0.5 Deep learning (splicing grammar) Jaganathan et al. 2019 Cell
SQUIRLS Splicing 0–1 ≥0.5 Graph-based splicing Riepe et al. 2021 AJHG
EVE Missense 0–1 ≥0.5 VAE evolutionary model Frazer et al. 2021 Nature
PrimateAI Missense 0–1 ≥0.803 Deep CNN on primate variation Sundaram et al. 2018 Nature Genetics

6.3 ACMG/AMP Variant Classification Framework

The ACMG/AMP framework (Richards et al. 2015) is the international standard for classifying germline variants in clinical contexts. It uses a two-letter code system with evidence strength (Very Strong, Strong, Moderate, Supporting):

Pathogenic evidence categories:
- PVS1 — Null variant (LOF) in gene where LOF is a disease mechanism (Very Strong)
- PS1 — Same AA change as established pathogenic variant (Strong)
- PM2 — Absent/very rare in gnomAD (Moderate)
- PP3 — Multiple computational predictors suggest damaging (Supporting)

Benign evidence categories:
- BA1 — Allele frequency >5% in gnomAD (Standalone Benign)
- BS1 — Allele frequency >expected for disorder (Strong Benign)
- BP4 — Multiple computational predictors suggest benign (Supporting Benign)

Reference: Richards S. et al. (2015). “Standards and guidelines for the interpretation of sequence variants.” Genetics in Medicine 17:405–424. [Link]

# R: implement basic ACMG classification score
classify_acmg <- function(pvs=0, ps=0, pm=0, pp=0, ba=0, bs=0, bp=0) {
  # Points-based approximation (Tavtigian et al. 2020 modification)
  path_score <- pvs*8 + ps*4 + pm*2 + pp*1
  benign_score <- ba*8 + bs*4 + bp*1
  
  net <- path_score - benign_score
  
  classification <- case_when(
    ba >= 1                    ~ "Benign (BA1 standalone)",
    net >= 18                  ~ "Pathogenic",
    net >= 10                  ~ "Likely Pathogenic",
    net >= 6  & net < 10       ~ "Variant of Uncertain Significance (VUS) — leaning Pathogenic",
    net >= -5 & net <= 5       ~ "Variant of Uncertain Significance (VUS)",
    net <= -7                  ~ "Likely Benign",
    net <= -8                  ~ "Benign",
    TRUE                       ~ "Variant of Uncertain Significance (VUS)"
  )
  
  cat(sprintf("Pathogenic evidence score: %d\n", path_score))
  cat(sprintf("Benign evidence score:     %d\n", benign_score))
  cat(sprintf("Net score:                 %d\n", net))
  cat(sprintf("Classification:            %s\n", classification))
  return(classification)
}

# Example: BRCA1 frameshift (PVS1) + absent in gnomAD (PM2) + splicing predictor (PP3)
cat("=== Example: BRCA1 frameshift variant ===\n")
#> === Example: BRCA1 frameshift variant ===
classify_acmg(pvs=1, pm=1, pp=1)
#> Pathogenic evidence score: 11
#> Benign evidence score:     0
#> Net score:                 11
#> Classification:            Likely Pathogenic
#> [1] "Likely Pathogenic"
cat("\n=== Example: Common missense, benign predictors ===\n")
#> 
#> === Example: Common missense, benign predictors ===
classify_acmg(pm=1, ba=1)
#> Pathogenic evidence score: 2
#> Benign evidence score:     8
#> Net score:                 -6
#> Classification:            Benign (BA1 standalone)
#> [1] "Benign (BA1 standalone)"

6.4 Key Variant Databases

dbs <- data.frame(
  Database = c("gnomAD v4","ClinVar","OMIM","dbSNP","HGMD (subscription)","ClinGen",
               "LOVD","UK Biobank","TCGA","COSMIC v99","InterVar","Franklin (Genoox)"),
  Content = c("Population allele frequencies (800K+ genomes/exomes)",
               "Variant-disease classifications (NCI/NIH hosted)",
               "Mendelian disease-gene relationships",
               "All submitted variants (900M+); rsID assignment",
               "Disease-causing mutations (curated)",
               "Gene-disease validity curation",
               "Locus-specific variant databases",
               "500K UK population WGS + phenotypes",
               "Pan-cancer somatic mutations (33 cancer types)",
               "Somatic mutation catalogue (~150M variants)",
               "ACMG classification software",
               "AI-assisted clinical variant interpretation"),
  URL = c("gnomad.broadinstitute.org","clinvar.ncbi.nlm.nih.gov",
          "omim.org","ncbi.nlm.nih.gov/snp",
          "hgmd.cf.ac.uk","clinicalgenome.org",
          "lovd.nl","ukbiobank.ac.uk",
          "portal.gdc.cancer.gov","cancer.sanger.ac.uk/cosmic",
          "intervar.missionbio.com","franklin.genoox.com"),
  stringsAsFactors=FALSE
)
kable(dbs, caption="Key variant and clinical databases") %>%
  kable_styling(bootstrap_options=c("striped","hover","condensed"), full_width=TRUE, font_size=11) %>%
  column_spec(1, bold=TRUE)
Key variant and clinical databases
Database Content URL
gnomAD v4 Population allele frequencies (800K+ genomes/exomes) gnomad.broadinstitute.org
ClinVar Variant-disease classifications (NCI/NIH hosted) clinvar.ncbi.nlm.nih.gov
OMIM Mendelian disease-gene relationships omim.org
dbSNP All submitted variants (900M+); rsID assignment ncbi.nlm.nih.gov/snp
HGMD (subscription) Disease-causing mutations (curated) hgmd.cf.ac.uk
ClinGen Gene-disease validity curation clinicalgenome.org
LOVD Locus-specific variant databases lovd.nl
UK Biobank 500K UK population WGS + phenotypes ukbiobank.ac.uk
TCGA Pan-cancer somatic mutations (33 cancer types) portal.gdc.cancer.gov
COSMIC v99 Somatic mutation catalogue (~150M variants) cancer.sanger.ac.uk/cosmic
InterVar ACMG classification software intervar.missionbio.com
Franklin (Genoox) AI-assisted clinical variant interpretation franklin.genoox.com

7 PART 5 — Bulk RNA-seq: From Raw Reads to Biological Insight

7.1 Library Preparation Strategies

Poly-A selection vs. rRNA depletion:
- Poly-A selection (TruSeq): uses oligo-dT beads to capture only polyadenylated (mature mRNA) transcripts. Pros: clean signal, less depth required. Cons: misses non-coding RNAs, snRNA, enhancer RNA; degrades with low-quality RNA.
- rRNA depletion (RiboMinus/Ribo-Zero): removes abundant rRNA (>80% of total RNA) and retains everything else. Pros: captures lncRNA, circRNA, pre-mRNA, viral RNA. Cons: more depth required; background noisier.
- stranded libraries: preserve strand information (critical for antisense transcription analysis and overlapping gene loci) — use --strandedness RF in STAR.

7.2 Splice-Aware Alignment with STAR

# Step 1: Generate genome index (one-time, takes 30-60 min, 32GB RAM)
STAR \
  --runMode genomeGenerate \
  --genomeDir star_index/ \
  --genomeFastaFiles GRCh38.fa \
  --sjdbGTFfile gencode.v44.annotation.gtf \
  --sjdbOverhang 149 \        # read length - 1; critical for splice junction detection
  --runThreadN 16

# Step 2: Align (2-pass mode for better novel junction discovery)
STAR \
  --genomeDir star_index/ \
  --readFilesIn trim_R1.fastq.gz trim_R2.fastq.gz \
  --readFilesCommand zcat \
  --outSAMtype BAM SortedByCoordinate \
  --outSAMattributes NH HI AS NM MD \
  --outSAMstrandField intronMotif \
  --quantMode GeneCounts \
  --sjdbGTFfile gencode.v44.annotation.gtf \
  --twopassMode Basic \
  --runThreadN 16 \
  --limitBAMsortRAM 32000000000 \
  --outFileNamePrefix results/sample1/

# STAR output files:
# Aligned.sortedByCoord.out.bam  — alignments
# ReadsPerGene.out.tab            — gene-level counts (column 2: unstranded, 3: sense, 4: antisense)
# SJ.out.tab                      — splice junction table
# Log.final.out                   — alignment statistics (check %uniquely mapped reads!)

Industry Q (Genentech, AstraZeneca RNA-seq teams): “A sample shows only 60% uniquely mapped reads in STAR — what could cause this and how do you investigate?”

Answer: Expected: >80% uniquely mapped for good RNA-seq data. Causes of low unique mapping:
1. rRNA contamination — STAR aligns rRNA reads but marks them as multiple-mapper; check featureCounts rRNA fraction or use RSeQC infer_experiment.py
2. Species contamination — BLAST overrepresented unmapped sequences
3. Library prep failure — adapter dimers, genomic DNA contamination
4. Wrong genome/annotation — check genome version matches species/annotation version
5. Strand library mismatch — wrong --strandedness setting inflates multi-mappers

Diagnostic: run RSeQC bam_stat.py, examine STAR’s Log.final.out, and BLAST the sequences in Unmapped.out.mate1.fastq.

7.3 Salmon — Alignment-Free Quantification

# Build index on transcriptome (decoy-aware — improves accuracy ~5%)
grep "^>" GRCh38.fa | cut -d " " -f 1 | sed 's/>//' > decoy_list.txt
cat transcriptome.fa GRCh38.fa > gentrome.fa

salmon index \
  -t gentrome.fa \
  -d decoy_list.txt \
  -i salmon_index \
  -p 16

# Quantify
salmon quant \
  -i salmon_index \
  -l A \
  -1 trim_R1.fastq.gz -2 trim_R2.fastq.gz \
  -p 16 \
  --validateMappings \
  --gcBias \
  --seqBias \
  -o quant_output/

# Output: quant.sf — columns: Name, Length, EffectiveLength, TPM, NumReads
# R: import and visualize Salmon quantification output
# Simulated quant.sf data
set.seed(42)
n_genes <- 200
salmon_sim <- data.frame(
  Name = paste0("ENST", sprintf("%011d", 1:n_genes), ".1"),
  Gene = paste0("GENE", 1:n_genes),
  Length = sample(500:5000, n_genes, replace=TRUE),
  TPM_treated = c(rlnorm(180, 3, 2), rep(0, 20)),
  TPM_control = c(rlnorm(170, 3, 2), rep(0, 30)),
  stringsAsFactors=FALSE
)
salmon_sim$TPM_treated[1:10] <- salmon_sim$TPM_treated[1:10] * 5  # simulate DE genes

# TPM vs length (should show no correlation post-length-normalization)
p1 <- ggplot(salmon_sim %>% filter(TPM_treated > 0), 
             aes(x=log10(Length), y=log10(TPM_treated+0.01))) +
  geom_point(alpha=0.4, color="#3498db") +
  geom_smooth(method="lm", color="red", se=FALSE) +
  labs(title="TPM vs Transcript Length\n(post-normalization, expect flat)", 
       x="log10(Length)", y="log10(TPM)") +
  theme_minimal(base_size=11)

p2 <- ggplot(salmon_sim, aes(x=log10(TPM_treated+0.01), y=log10(TPM_control+0.01))) +
  geom_point(alpha=0.4, aes(color=TPM_treated/pmax(TPM_control,0.01) > 4)) +
  geom_abline(slope=1, intercept=0, color="red", linetype="dashed") +
  scale_color_manual(values=c("FALSE"="#bdc3c7","TRUE"="#e74c3c"),
                     labels=c("Not DE","Potential upregulation")) +
  labs(title="MA-style: Treated vs Control TPM", 
       x="log10(TPM treated)", y="log10(TPM control)", color="") +
  theme_minimal(base_size=11)

grid.arrange(p1, p2, ncol=2)

Reference: Patro R. et al. (2017). “Salmon provides fast and bias-aware quantification of transcript expression.” Nature Methods 14:417–419. [Link]

7.4 DESeq2 — Differential Expression Analysis

# Complete DESeq2 workflow with simulated RNA-seq data
# In production: replace with actual count matrix from featureCounts/STAR GeneCounts/tximeta

set.seed(123)
n_genes <- 1000
n_samples <- 6  # 3 treated, 3 control

# Simulate count data (negative binomial)
# Real pipeline: dds <- DESeqDataSetFromTximeta(tximeta_output, design = ~condition)
sim_counts <- matrix(
  c(rnbinom(n_genes * 3, mu=100, size=10),     # control samples
    rnbinom(n_genes * 3, mu=c(rep(100, n_genes*2), rep(500, n_genes)), size=10)),  # treated
  nrow=n_genes, ncol=n_samples
)
# Add DE signal to first 50 genes
sim_counts[1:50, 4:6] <- sim_counts[1:50, 4:6] * sample(c(3,4,5), 50, replace=TRUE)
sim_counts[51:100, 4:6] <- round(sim_counts[51:100, 4:6] * 0.2)  # downregulated

rownames(sim_counts) <- paste0("ENSG", sprintf("%011d", 1:n_genes))
colnames(sim_counts) <- c(paste0("ctrl_", 1:3), paste0("treat_", 1:3))

# Column data
coldata <- data.frame(
  condition = factor(c(rep("control", 3), rep("treated", 3))),
  batch = factor(c("A","A","B","A","B","B")),
  row.names = colnames(sim_counts)
)

cat("=== Simulated Count Matrix (first 5 genes, 6 samples) ===\n")
#> === Simulated Count Matrix (first 5 genes, 6 samples) ===
print(head(sim_counts, 5))
#>                 ctrl_1 ctrl_2 ctrl_3 treat_1 treat_2 treat_3
#> ENSG00000000001     89    164    102     395     480    1330
#> ENSG00000000002     50     77     98     387     258    2205
#> ENSG00000000003    160     57     87     560     270    3115
#> ENSG00000000004     46     97    145     436     384    1860
#> ENSG00000000005    140    107    101    1005     380    3760
cat("\n=== Sample metadata ===\n")
#> 
#> === Sample metadata ===
print(coldata)
#>         condition batch
#> ctrl_1    control     A
#> ctrl_2    control     A
#> ctrl_3    control     B
#> treat_1   treated     A
#> treat_2   treated     B
#> treat_3   treated     B
if (requireNamespace("DESeq2", quietly=TRUE)) {
  library(DESeq2)
  
  # Build DESeq2 object
  dds <- DESeqDataSetFromMatrix(
    countData = sim_counts,
    colData   = coldata,
    design    = ~ batch + condition   # account for batch effect!
  )
  
  # Pre-filter: remove genes with <10 counts total
  keep <- rowSums(counts(dds)) >= 10
  dds  <- dds[keep, ]
  cat(sprintf("Genes retained after low-count filter: %d / %d\n", sum(keep), nrow(sim_counts)))
  
  # Run DESeq2 (size-factor estimation + dispersion modelling + Wald test)
  dds <- DESeq(dds)
  
  # Extract results with log2FC shrinkage (apeglm — recommended)
  res <- lfcShrink(dds, coef="condition_treated_vs_control", type="apeglm")
  res_df <- as.data.frame(res)
  res_df$gene <- rownames(res_df)
  res_df$significant <- !is.na(res_df$padj) & res_df$padj < 0.05
  res_df$direction <- ifelse(res_df$log2FoldChange > 0, "Up", "Down")
  
  cat("\n=== DESeq2 Results Summary ===\n")
  print(summary(res, alpha=0.05))
  
  # Volcano plot
  volcano_df <- res_df %>% filter(!is.na(padj), !is.na(log2FoldChange))
  
  p_volcano <- ggplot(volcano_df, aes(x=log2FoldChange, y=-log10(pvalue), 
                                       color=significant & abs(log2FoldChange) > 1)) +
    geom_point(alpha=0.5, size=1.5) +
    geom_vline(xintercept=c(-1,1), linetype="dashed", color="grey40") +
    geom_hline(yintercept=-log10(0.05), linetype="dashed", color="grey40") +
    scale_color_manual(values=c("FALSE"="#bdc3c7","TRUE"="#e74c3c"),
                       labels=c("FALSE"="Not significant","TRUE"="DE (padj<0.05, |LFC|>1)")) +
    geom_text_repel(data=volcano_df %>% filter(significant, abs(log2FoldChange)>2) %>% head(10),
                    aes(label=gene), size=2.5, max.overlaps=10) +
    labs(title="Volcano Plot: Treated vs Control",
         subtitle="DESeq2 with apeglm LFC shrinkage",
         x="Shrunken log2 Fold Change", y="-log10(p-value)", color="") +
    theme_minimal(base_size=12) + theme(legend.position="bottom")
  
  # MA plot
  p_ma <- ggplot(volcano_df, aes(x=log10(baseMean+1), y=log2FoldChange,
                                   color=significant)) +
    geom_point(alpha=0.4, size=1.2) +
    geom_hline(yintercept=0, color="black") +
    scale_color_manual(values=c("FALSE"="#bdc3c7","TRUE"="#e74c3c"), guide=FALSE) +
    labs(title="MA Plot", x="log10(Mean Expression)", y="log2 Fold Change") +
    theme_minimal(base_size=12)
  
  grid.arrange(p_volcano, p_ma, ncol=2)
  
  # Top DE genes table
  top_de <- res_df %>%
    filter(significant) %>%
    arrange(padj) %>%
    head(15) %>%
    select(gene, baseMean, log2FoldChange, lfcSE, pvalue, padj) %>%
    mutate(across(where(is.numeric), ~round(., 4)))
  
  kable(top_de, caption="Top 15 Differentially Expressed Genes (DESeq2, padj<0.05)") %>%
    kable_styling(bootstrap_options=c("striped","hover"), full_width=FALSE)
  
} else {
  cat("Install DESeq2: BiocManager::install('DESeq2')\n")
  # Show mock results table
  mock_results <- data.frame(
    gene = paste0("ENSG", sprintf("%011d", 1:10)),
    baseMean = round(runif(10, 50, 2000), 1),
    log2FoldChange = round(c(3.2, 2.8, -2.1, 2.4, -3.1, 1.9, -1.8, 2.6, -2.3, 3.0), 2),
    padj = signif(c(1e-15, 2e-12, 5e-10, 3e-9, 1e-8, 4e-7, 2e-6, 5e-6, 8e-6, 1e-5), 2)
  )
  kable(mock_results, caption="Example DESeq2 results (mock data)") %>%
    kable_styling(bootstrap_options="striped", full_width=FALSE)
}
#> Genes retained after low-count filter: 1000 / 1000
#> 
#> === DESeq2 Results Summary ===
#> 
#> out of 1000 with nonzero total read count
#> adjusted p-value < 0.05
#> LFC > 0 (up)       : 53, 5.3%
#> LFC < 0 (down)     : 52, 5.2%
#> outliers [1]       : 0, 0%
#> low counts [2]     : 0, 0%
#> (mean count < 54)
#> [1] see 'cooksCutoff' argument of ?results
#> [2] see 'independentFiltering' argument of ?results
#> 
#> NULL
Top 15 Differentially Expressed Genes (DESeq2, padj<0.05)
gene baseMean log2FoldChange lfcSE pvalue padj
ENSG00000000018 ENSG00000000018 428.1865 2.7365 0.3530 0 0
ENSG00000000043 ENSG00000000043 380.0211 2.8834 0.3636 0 0
ENSG00000000024 ENSG00000000024 418.3991 2.7222 0.3679 0 0
ENSG00000000094 ENSG00000000094 96.7858 -3.4563 0.4769 0 0
ENSG00000000084 ENSG00000000084 96.1694 -2.8566 0.4086 0 0
ENSG00000000025 ENSG00000000025 411.7784 2.5735 0.3681 0 0
ENSG00000000005 ENSG00000000005 531.8095 2.6374 0.3954 0 0
ENSG00000000054 ENSG00000000054 79.4543 -2.7841 0.4552 0 0
ENSG00000000068 ENSG00000000068 82.4993 -2.6455 0.4201 0 0
ENSG00000000077 ENSG00000000077 84.8032 -2.7389 0.4301 0 0
ENSG00000000091 ENSG00000000091 81.6750 -2.7102 0.4294 0 0
ENSG00000000048 ENSG00000000048 376.2193 2.2789 0.3562 0 0
ENSG00000000034 ENSG00000000034 381.4179 2.4110 0.4227 0 0
ENSG00000000096 ENSG00000000096 90.9570 -2.5446 0.4453 0 0
ENSG00000000087 ENSG00000000087 93.3655 -2.4037 0.3973 0 0

References: Love M.I., Huber W., Anders S. (2014). “Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2.” Genome Biology 15:550. [Link]
Zhu A., Ibrahim J.G., Love M.I. (2019). “Heavy-tailed prior distributions for sequence count data: removing the noise and preserving large differences.” Bioinformatics 35(12):2084–2092. [apeglm Link]

7.5 Gene Set Enrichment Analysis (GSEA) and ORA

if (requireNamespace("clusterProfiler", quietly=TRUE) && 
    requireNamespace("org.Hs.eg.db", quietly=TRUE)) {
  
  library(clusterProfiler)
  library(org.Hs.eg.db)
  
  # Gene Ontology Over-Representation Analysis (ORA)
  # In real pipeline: use significant DE genes
  sig_genes <- paste0("ENSG", sprintf("%011d", sample(1:20000, 200)))
  
  # Convert to Entrez IDs
  gene_map <- bitr(sig_genes, fromType="ENSEMBL", toType="ENTREZID", OrgDb=org.Hs.eg.db)
  
  ego <- enrichGO(
    gene          = gene_map$ENTREZID,
    OrgDb         = org.Hs.eg.db,
    ont           = "BP",  # Biological Process
    pAdjustMethod = "BH",
    qvalueCutoff  = 0.05,
    readable      = TRUE
  )
  
  if (!is.null(ego) && nrow(as.data.frame(ego)) > 0) {
    dotplot(ego, showCategory=15) + 
      labs(title="GO Biological Process Enrichment (ORA)") +
      theme_minimal(base_size=11)
  }
} else {
  # Visualize simulated enrichment results
  enrich_sim <- data.frame(
    GO_Term = c("DNA repair","Cell cycle checkpoint","mRNA processing",
                 "Apoptosis regulation","Immune response","Metabolic process",
                 "Transcription regulation","Signal transduction","Cell migration",
                 "Chromatin remodeling"),
    neg_log10_padj = c(12.5, 10.2, 8.8, 7.3, 6.9, 5.4, 4.8, 4.2, 3.7, 3.1),
    Gene_Ratio = c(0.35, 0.28, 0.31, 0.22, 0.19, 0.41, 0.27, 0.33, 0.18, 0.24),
    Count = c(42, 35, 38, 28, 24, 51, 34, 41, 22, 30)
  )
  
  ggplot(enrich_sim, aes(x=Gene_Ratio, y=reorder(GO_Term, neg_log10_padj),
                          color=neg_log10_padj, size=Count)) +
    geom_point() +
    scale_color_gradient(low="#3498db", high="#e74c3c") +
    labs(title="GO Biological Process Enrichment (Simulated)",
         x="Gene Ratio", y="GO Term",
         color="-log10(padj)", size="Gene Count") +
    theme_minimal(base_size=12)
}

Reference: Subramanian A. et al. (2005). “Gene set enrichment analysis.” PNAS 102(43):15545–15550. [Link]
Yu G. et al. (2012). “clusterProfiler: an R Package for Comparing Biological Themes Among Gene Clusters.” OMICS 16(5):284–287. [Link]

8 PART 6 — Single-Cell RNA-seq: Complete Workflow

8.1 10x Genomics Chemistry — Read Architecture

# Visualise 10x v3 read structure
read_parts <- data.frame(
  Component = c("Cell Barcode (16 bp)","UMI (12 bp)","Poly-dT / TSO","cDNA insert (~90 bp)"),
  Read = c("R1","R1","R1","R2"),
  Start_bp = c(1, 17, 29, 1),
  Width_bp = c(16, 12, 12, 90),
  Color = c("#e74c3c","#3498db","#f39c12","#27ae60"),
  stringsAsFactors=FALSE
)

ggplot() +
  # R1
  annotate("text", x=-2, y=2, label="R1 (28bp)", hjust=1, fontface="bold", size=4) +
  geom_rect(data=read_parts %>% filter(Read=="R1"),
            aes(xmin=Start_bp, xmax=Start_bp+Width_bp, ymin=1.6, ymax=2.4, fill=Component)) +
  geom_text(data=read_parts %>% filter(Read=="R1"),
            aes(x=Start_bp+Width_bp/2, y=2, label=paste0(Component,"\n",Width_bp,"bp")),
            size=3, fontface="bold", color="white") +
  # R2
  annotate("text", x=-2, y=1, label="R2 (~90bp)", hjust=1, fontface="bold", size=4) +
  geom_rect(data=read_parts %>% filter(Read=="R2"),
            aes(xmin=Start_bp, xmax=Start_bp+Width_bp, ymin=0.6, ymax=1.4, fill=Component)) +
  geom_text(data=read_parts %>% filter(Read=="R2"),
            aes(x=Start_bp+Width_bp/2, y=1, label=paste0(Component,"\n",Width_bp,"bp")),
            size=3, fontface="bold", color="white") +
  scale_fill_manual(values=setNames(read_parts$Color, read_parts$Component)) +
  scale_x_continuous(breaks=seq(0,110,10)) +
  coord_cartesian(xlim=c(-5, 110), ylim=c(0.2, 2.8)) +
  labs(title="10x Genomics Chromium v3 Read Structure",
       subtitle="R1 carries the cell barcode + UMI; R2 carries the cDNA sequence",
       x="Position in read (bp)", y="", fill="Component") +
  theme_minimal(base_size=12) +
  theme(axis.text.y=element_blank(), axis.ticks.y=element_blank())

8.2 STARsolo / Cell Ranger — Count Matrix Generation

# STARsolo (open-source alternative to Cell Ranger, 3x faster)
STAR \
  --soloType CB_UMI_Simple \
  --soloCBwhitelist 3M-february-2018.txt \
  --soloCBstart 1 --soloCBlen 16 \
  --soloUMIstart 17 --soloUMIlen 12 \
  --readFilesIn R2.fastq.gz R1.fastq.gz \     # Note: cDNA first, barcode second!
  --readFilesCommand zcat \
  --genomeDir star_index/ \
  --outSAMtype BAM SortedByCoordinate \
  --outSAMattributes NH HI nM AS CR UR CB UB GX GN sS sQ sM \
  --soloFeatures Gene SJ \
  --soloCellFilter EmptyDrops_CR \
  --soloOutDir Solo.out/ \
  --runThreadN 16

# Cell Ranger (10x proprietary — wrapper around STAR + custom algorithms)
cellranger count \
  --id=sample1 \
  --transcriptome=/path/to/refdata-gex-GRCh38-2024-A \
  --fastqs=/path/to/fastqs \
  --sample=sample1 \
  --localcores=32 \
  --localmem=128 \
  --expect-cells=5000

References: Kaminow B., Yunusov D., Dobin A. (2021). STARsolo. bioRxiv 2021.05.05.442755; 10x Genomics Cell Ranger documentation. [Cell Ranger]

8.3 Complete Scanpy/Seurat Analysis Pipeline

# ====================================================================
# Complete scRNA-seq analysis pipeline in Python (Scanpy)
# ====================================================================
import scanpy as sc
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')

sc.settings.verbosity = 1
sc.settings.figdir = "figures/"

# ── 1. Load data ────────────────────────────────────────────────────
adata = sc.read_10x_h5("filtered_feature_bc_matrix.h5")
adata.var_names_make_unique()
print(f"Loaded: {adata.n_obs} cells × {adata.n_vars} genes")

# ── 2. Quality metrics ──────────────────────────────────────────────
adata.var['mt'] = adata.var_names.str.startswith('MT-')
adata.var['ribo'] = adata.var_names.str.startswith(('RPS','RPL'))
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt','ribo'], 
                            percent_top=None, log1p=False, inplace=True)

print("\nQC summary:")
print(f"  Median genes/cell:    {adata.obs.n_genes_by_counts.median():.0f}")
print(f"  Median UMIs/cell:     {adata.obs.total_counts.median():.0f}")
print(f"  Median %MT:           {adata.obs.pct_counts_mt.median():.1f}%")

# Visualise QC
fig, axes = plt.subplots(1, 3, figsize=(12, 4))
for ax, metric, color in zip(axes,
    ['n_genes_by_counts','total_counts','pct_counts_mt'],
    ['steelblue','salmon','lightgreen']):
    ax.hist(adata.obs[metric], bins=100, color=color, edgecolor='none')
    ax.set_xlabel(metric.replace('_',' '))
    ax.set_ylabel('# cells')
    ax.set_title(metric.replace('_',' ').title())
plt.tight_layout()
plt.savefig("figures/qc_histograms.png", dpi=150, bbox_inches='tight')

# ── 3. Filtering ────────────────────────────────────────────────────
MIN_GENES = 200;  MAX_GENES = 6000   # remove doublets (high) and empty drops (low)
MAX_MT = 15;      MIN_COUNTS = 500
adata = adata[
    (adata.obs.n_genes_by_counts > MIN_GENES) &
    (adata.obs.n_genes_by_counts < MAX_GENES) &
    (adata.obs.pct_counts_mt < MAX_MT) &
    (adata.obs.total_counts > MIN_COUNTS)
].copy()
sc.pp.filter_genes(adata, min_cells=10)  # remove genes in <10 cells
print(f"\nAfter QC filtering: {adata.n_obs} cells × {adata.n_vars} genes")

# ── 4. Doublet detection with Scrublet ──────────────────────────────
try:
    import scrublet as scr
    scrub = scr.Scrublet(adata.X, expected_doublet_rate=0.08)
    doublet_scores, predicted_doublets = scrub.scrub_doublets(min_counts=3, min_cells=3)
    adata.obs['doublet_score'] = doublet_scores
    adata.obs['predicted_doublet'] = predicted_doublets
    adata = adata[~adata.obs.predicted_doublet].copy()
    print(f"After doublet removal: {adata.n_obs} cells")
except ImportError:
    print("Scrublet not installed: pip install scrublet")

# ── 5. Normalisation ────────────────────────────────────────────────
adata.layers['counts'] = adata.X.copy()  # keep raw counts for DE

# Scran-pooling normalisation (gold standard, requires rpy2 + scran R package)
# Simplified: use log-normalisation as fallback
sc.pp.normalize_total(adata, target_sum=1e4)    # scale each cell to 10,000 counts
sc.pp.log1p(adata)                               # log(X + 1)
adata.raw = adata                                # freeze normalised counts for DE

# ── 6. Feature selection ────────────────────────────────────────────
sc.pp.highly_variable_genes(adata, n_top_genes=3000, subset=True, flavor='seurat_v3')
print(f"Highly variable genes selected: {adata.n_vars}")

# ── 7. Scaling (for PCA) ────────────────────────────────────────────
sc.pp.scale(adata, max_value=10)  # Z-score per gene, clip at 10 SD

# ── 8. PCA ──────────────────────────────────────────────────────────
sc.tl.pca(adata, svd_solver='arpack', n_comps=50)
sc.pl.pca_variance_ratio(adata, n_pcs=50, log=True, save='_variance_ratio.png')

# ── 9. Batch integration (Harmony) ──────────────────────────────────
try:
    sc.external.pp.harmony_integrate(adata, key='batch', basis='X_pca', 
                                      adjusted_basis='X_pca_harmony')
    use_rep = 'X_pca_harmony'
    print("Harmony integration applied")
except Exception:
    use_rep = 'X_pca'
    print("Harmony not available — using PCA directly")

# ── 10. Neighbourhood graph, UMAP, clustering ───────────────────────
sc.pp.neighbors(adata, n_neighbors=20, n_pcs=30, use_rep=use_rep)
sc.tl.umap(adata, min_dist=0.3, spread=1.0)
sc.tl.leiden(adata, resolution=0.6, random_state=42)

print(f"\nLeiden clustering: {adata.obs.leiden.nunique()} clusters")

# ── 11. Marker gene detection ───────────────────────────────────────
sc.tl.rank_genes_groups(adata, 'leiden', method='wilcoxon', 
                          pts=True, use_raw=True)

# ── 12. Canonical cell-type markers ─────────────────────────────────
canonical_markers = {
    'T cells (Pan)': ['CD3D','CD3E','CD3G'],
    'CD8+ T cells': ['CD8A','CD8B'],
    'CD4+ T cells': ['CD4','IL7R'],
    'Tregs': ['FOXP3','IL2RA','CTLA4'],
    'NK cells': ['NKG7','GNLY','KLRD1'],
    'B cells': ['MS4A1','CD79A','CD19'],
    'Plasma cells': ['MZB1','JCHAIN','SDC1'],
    'TAM (M2-like)': ['CD68','CD163','CSF1R'],
    'CAFs': ['COL1A1','PDGFRB','FAP'],
    'Endothelial': ['PECAM1','VWF','CDH5'],
    'Epithelial/Tumor': ['EPCAM','KRT18','KRT8']
}

# Compute module scores for each lineage
for cell_type, markers in canonical_markers.items():
    valid_markers = [g for g in markers if g in adata.var_names]
    if valid_markers:
        sc.tl.score_genes(adata, valid_markers, score_name=f"score_{cell_type.split()[0]}")

print("\nCell type scoring complete")
print(adata.obs[[c for c in adata.obs.columns if c.startswith('score_')]].describe().round(3))
# R equivalent — Seurat v5 pipeline
if (requireNamespace("Seurat", quietly=TRUE)) {
  library(Seurat)
  
  cat("=== Seurat v5 pipeline outline ===\n")
  cat("
# Load 10x data
counts <- Read10X(data.dir = 'filtered_feature_bc_matrix/')
sobj   <- CreateSeuratObject(counts=counts, min.cells=3, min.features=200, 
                               project='TumorSample')

# QC
sobj[['percent.mt']]   <- PercentageFeatureSet(sobj, pattern='^MT-')
sobj[['percent.ribo']] <- PercentageFeatureSet(sobj, pattern='^RP[SL]')

# Filter
sobj <- subset(sobj, 
               nFeature_RNA > 200 & nFeature_RNA < 6000 & 
               percent.mt < 15 & nCount_RNA > 500)

# Normalise + Highly variable genes + Scale (SCTransform replaces all 3)
sobj <- SCTransform(sobj, vars.to.regress=c('percent.mt','nCount_RNA'),
                    variable.features.n=3000, verbose=FALSE)

# Dimensionality reduction
sobj <- RunPCA(sobj, npcs=50)

# Integration (Harmony)
library(harmony)
sobj <- RunHarmony(sobj, group.by.vars='batch', reduction.use='pca',
                   reduction.save='harmony')

# Clustering + UMAP
sobj <- FindNeighbors(sobj, reduction='harmony', dims=1:30)
sobj <- FindClusters(sobj, resolution=0.6, algorithm=1)  # Louvain
sobj <- RunUMAP(sobj, reduction='harmony', dims=1:30)

# Marker genes
markers <- FindAllMarkers(sobj, only.pos=TRUE, min.pct=0.25,
                           logfc.threshold=0.25, test.use='wilcox')

# Top markers per cluster
top_markers <- markers %>% group_by(cluster) %>% slice_max(avg_log2FC, n=5)
print(top_markers)
")
} else {
  cat("Install Seurat: install.packages('Seurat')\n")
}
#> === Seurat v5 pipeline outline ===
#> 
#> # Load 10x data
#> counts <- Read10X(data.dir = 'filtered_feature_bc_matrix/')
#> sobj   <- CreateSeuratObject(counts=counts, min.cells=3, min.features=200, 
#>                                project='TumorSample')
#> 
#> # QC
#> sobj[['percent.mt']]   <- PercentageFeatureSet(sobj, pattern='^MT-')
#> sobj[['percent.ribo']] <- PercentageFeatureSet(sobj, pattern='^RP[SL]')
#> 
#> # Filter
#> sobj <- subset(sobj, 
#>                nFeature_RNA > 200 & nFeature_RNA < 6000 & 
#>                percent.mt < 15 & nCount_RNA > 500)
#> 
#> # Normalise + Highly variable genes + Scale (SCTransform replaces all 3)
#> sobj <- SCTransform(sobj, vars.to.regress=c('percent.mt','nCount_RNA'),
#>                     variable.features.n=3000, verbose=FALSE)
#> 
#> # Dimensionality reduction
#> sobj <- RunPCA(sobj, npcs=50)
#> 
#> # Integration (Harmony)
#> library(harmony)
#> sobj <- RunHarmony(sobj, group.by.vars='batch', reduction.use='pca',
#>                    reduction.save='harmony')
#> 
#> # Clustering + UMAP
#> sobj <- FindNeighbors(sobj, reduction='harmony', dims=1:30)
#> sobj <- FindClusters(sobj, resolution=0.6, algorithm=1)  # Louvain
#> sobj <- RunUMAP(sobj, reduction='harmony', dims=1:30)
#> 
#> # Marker genes
#> markers <- FindAllMarkers(sobj, only.pos=TRUE, min.pct=0.25,
#>                            logfc.threshold=0.25, test.use='wilcox')
#> 
#> # Top markers per cluster
#> top_markers <- markers %>% group_by(cluster) %>% slice_max(avg_log2FC, n=5)
#> print(top_markers)

8.4 scRNA-seq QC Thresholds and Metrics

scrna_qc <- data.frame(
  Metric = c("Genes detected per cell","UMIs (counts) per cell","% Mitochondrial reads",
             "% Ribosomal reads","Doublet score (Scrublet)","Cells per library",
             "Median genes (10x v3)","Sequencing saturation"),
  Typical_Range = c("200–5,000 (tissue-dependent)","500–50,000","3–20%","5–40%","<0.2",
                     "1,000–15,000","~2,500–4,000 (PBMC)","60–80%"),
  Low_Warning = c("<200 (empty drops)","<500","—","—","—","<500 (poor capture)","<1,500","<50%"),
  High_Warning = c(">8,000 (doublets?)","High depends on cell type",">25% (dead cells)",">60%",
                   ">0.25 (high doublet risk)","—","—","—"),
  stringsAsFactors=FALSE
)
kable(scrna_qc, caption="scRNA-seq QC metrics and thresholds") %>%
  kable_styling(bootstrap_options=c("striped","hover","condensed"), full_width=TRUE) %>%
  column_spec(1, bold=TRUE)
scRNA-seq QC metrics and thresholds
Metric Typical_Range Low_Warning High_Warning
Genes detected per cell 200–5,000 (tissue-dependent) <200 (empty drops) >8,000 (doublets?)
UMIs (counts) per cell 500–50,000 <500 High depends on cell type
% Mitochondrial reads 3–20% — >25% (dead cells)
% Ribosomal reads 5–40% — >60%
Doublet score (Scrublet) <0.2 — >0.25 (high doublet risk)
Cells per library 1,000–15,000 <500 (poor capture) —
Median genes (10x v3) ~2,500–4,000 (PBMC) <1,500 —
Sequencing saturation 60–80% <50% —

9 PART 7 — Tumor Microenvironment (TME) Case Study

9.1 Expected Cell Types in Solid Tumors

tme_markers <- data.frame(
  Cell_Type = c("CD8+ Cytotoxic T","CD4+ Helper T","Tregs","NK cells","B cells",
                 "Plasma cells","M1-like TAM","M2-like TAM","cDC1","cDC2",
                 "CAFs","Endothelial","Malignant Epithelial","Mast cells","Neutrophils"),
  Lineage = c(rep("Lymphoid",5),"Lymphoid",rep("Myeloid",4),rep("Stromal",2),"Tumor","Myeloid","Myeloid"),
  Marker_Genes = c(
    "CD8A, CD8B, GZMB, PRF1, IFNG",
    "CD4, IL7R, CXCR4, MAL",
    "FOXP3, IL2RA, CTLA4, IKZF2",
    "NKG7, GNLY, KLRD1, NCAM1 (CD56; no CD3)",
    "MS4A1 (CD20), CD79A, CD19",
    "MZB1, JCHAIN, SDC1 (CD138), IGHG1",
    "CD68, CD80, CD86, CXCL10, HLA-DRA",
    "CD68, CD163, MRC1, ARG1, TGFB1",
    "CLEC9A, BATF3, XCR1, IRF8",
    "CLEC10A, LAMP3, CD1C, FCER1A",
    "COL1A1, COL1A2, PDGFRB, FAP, ACTA2",
    "PECAM1 (CD31), VWF, CDH5, CLDN5",
    "EPCAM, KRT18, KRT8, KRT19 + CNV inference",
    "KIT (CD117), TPSAB1, TPSB2, CPA3",
    "S100A8, S100A9, FCGR3B, CXCR2"
  ),
  Function_in_TME = c(
    "Anti-tumor cytotoxicity; exhaustion target of PD-1/TIM-3/LAG3",
    "Cytokine help; promote Treg or effector response",
    "Immunosuppression; FOXP3+ immunotherapy resistance marker",
    "Innate cytotoxicity; ADCC",
    "Antigen presentation; tertiary lymphoid structure formation",
    "Antibody secretion; often correlates with better prognosis",
    "Pro-inflammatory; promotes anti-tumor immunity",
    "Immunosuppression, angiogenesis; associated with poor outcome",
    "Cross-presentation; CD8 T cell priming (key for ICB response)",
    "CD4 T cell priming; bridge innate-adaptive",
    "ECM remodelling; drug delivery barrier; diverse states",
    "Tumour angiogenesis; immune cell trafficking; normalisation target",
    "Tumour cells; identified by aneuploidy via inferCNV/CopyKAT",
    "Tumour promotion; angiogenesis; recruitment target",
    "Pro- or anti-tumour; emerging immunotherapy target"
  ),
  stringsAsFactors=FALSE
)
kable(tme_markers, caption="Cell types, markers, and functions in the tumor microenvironment") %>%
  kable_styling(bootstrap_options=c("striped","hover","condensed"), full_width=TRUE, font_size=11) %>%
  column_spec(1, bold=TRUE) %>%
  column_spec(2, color="white", 
              background=ifelse(tme_markers$Lineage=="Lymphoid","#2980b9",
                               ifelse(tme_markers$Lineage=="Myeloid","#e67e22",
                                     ifelse(tme_markers$Lineage=="Tumor","#e74c3c","#27ae60"))))
Cell types, markers, and functions in the tumor microenvironment
Cell_Type Lineage Marker_Genes Function_in_TME
CD8+ Cytotoxic T Lymphoid CD8A, CD8B, GZMB, PRF1, IFNG Anti-tumor cytotoxicity; exhaustion target of PD-1/TIM-3/LAG3
CD4+ Helper T Lymphoid CD4, IL7R, CXCR4, MAL Cytokine help; promote Treg or effector response
Tregs Lymphoid FOXP3, IL2RA, CTLA4, IKZF2 Immunosuppression; FOXP3+ immunotherapy resistance marker
NK cells Lymphoid NKG7, GNLY, KLRD1, NCAM1 (CD56; no CD3) Innate cytotoxicity; ADCC
B cells Lymphoid MS4A1 (CD20), CD79A, CD19 Antigen presentation; tertiary lymphoid structure formation
Plasma cells Lymphoid MZB1, JCHAIN, SDC1 (CD138), IGHG1 Antibody secretion; often correlates with better prognosis
M1-like TAM Myeloid CD68, CD80, CD86, CXCL10, HLA-DRA Pro-inflammatory; promotes anti-tumor immunity
M2-like TAM Myeloid CD68, CD163, MRC1, ARG1, TGFB1 Immunosuppression, angiogenesis; associated with poor outcome
cDC1 Myeloid CLEC9A, BATF3, XCR1, IRF8 Cross-presentation; CD8 T cell priming (key for ICB response)
cDC2 Myeloid CLEC10A, LAMP3, CD1C, FCER1A CD4 T cell priming; bridge innate-adaptive
CAFs Stromal COL1A1, COL1A2, PDGFRB, FAP, ACTA2 ECM remodelling; drug delivery barrier; diverse states
Endothelial Stromal PECAM1 (CD31), VWF, CDH5, CLDN5 Tumour angiogenesis; immune cell trafficking; normalisation target
Malignant Epithelial Tumor EPCAM, KRT18, KRT8, KRT19 + CNV inference Tumour cells; identified by aneuploidy via inferCNV/CopyKAT
Mast cells Myeloid KIT (CD117), TPSAB1, TPSB2, CPA3 Tumour promotion; angiogenesis; recruitment target
Neutrophils Myeloid S100A8, S100A9, FCGR3B, CXCR2 Pro- or anti-tumour; emerging immunotherapy target

9.2 T Cell Exhaustion Scoring

# T cell exhaustion gene module scoring
exhaustion_genes <- c("PDCD1","HAVCR2","LAG3","CTLA4","TIGIT",  # inhibitory receptors
                        "ENTPD1","CD244","CD160","CXCL13")         # exhaustion markers

effector_genes   <- c("IFNG","GZMB","GZMH","PRF1","FASLG","TNF")  # effector function

# Simulate pseudotime + exhaustion/effector scores
set.seed(42)
n_cells <- 500
pseudotime <- seq(0, 1, length.out=n_cells) + rnorm(n_cells, 0, 0.05)
exhaustion_score <- pseudotime * 2 + rnorm(n_cells, 0, 0.3)
effector_score   <- (1 - pseudotime) * 2 + rnorm(n_cells, 0, 0.3)

state_labels <- cut(pseudotime,
                    breaks=c(-Inf, 0.3, 0.6, Inf),
                    labels=c("Effector/Progenitor","Intermediate","Terminally Exhausted"))

exh_df <- data.frame(Pseudotime=pseudotime, Exhaustion=exhaustion_score,
                      Effector=effector_score, State=state_labels)

p1 <- ggplot(exh_df, aes(x=Pseudotime, y=Exhaustion, color=State)) +
  geom_point(alpha=0.5, size=1) +
  geom_smooth(se=FALSE, method="loess", color="black", linewidth=1) +
  scale_color_brewer(palette="RdYlGn", direction=-1) +
  labs(title="Exhaustion Score Along CD8 T Cell Pseudotime",
       subtitle="PDCD1, HAVCR2 (TIM3), LAG3, CTLA4, TIGIT module",
       x="Pseudotime (naive → terminal exhaustion)", y="Exhaustion Score") +
  theme_minimal(base_size=11) + theme(legend.position="bottom")

p2 <- ggplot(exh_df, aes(x=Exhaustion, y=Effector, color=State)) +
  geom_point(alpha=0.5, size=1) +
  scale_color_brewer(palette="RdYlGn", direction=-1) +
  labs(title="Exhaustion vs. Effector Function Score",
       x="Exhaustion Score", y="Effector Score (GZMB, IFNG, PRF1)") +
  theme_minimal(base_size=11) + theme(legend.position="bottom")

grid.arrange(p1, p2, ncol=2)

References: Wherry E.J., Kurachi M. (2015). Nature Reviews Immunology 15:486–499; Miller B.C. et al. (2019). Nature Immunology 20:326–336.

10 PART 8 — Long-Read Sequencing: Nanopore and PacBio

10.1 Long-Read Applications

lr_apps <- data.frame(
  Application = c("Structural variant calling","Full-length isoform sequencing (Iso-Seq)",
                   "Direct RNA sequencing","De novo genome assembly",
                   "Repeat expansion genotyping","CpG methylation (direct)",
                   "Haplotype phasing","Metagenomics (species ID)"),
  Short_Read = c("Limited (Manta: >50bp only)","No (splice graphs only)",
                  "No (cDNA proxy)","No (fragmented assemblies)",
                  "Partial (ExpansionHunter)","No (bisulfite proxy only)",
                  "Limited (<10kb haplotypes)","Species-level (16S)"),
  PacBio_HiFi = c("Excellent (PBSV, Sniffles2)","Yes (Iso-Seq 3)",
                   "No","Best accuracy (hifiasm)","Yes (TRGT)",
                   "5mC with pb-CpG-tools","Yes, chromosome-scale","Yes (HiFi metagenomics)"),
  Oxford_Nanopore = c("Excellent (Sniffles2)","Yes (IsoSeq-like)","Yes (direct RNA)",
                       "Good (Flye, Verkko)","Yes (Straglr)","6mA, 5mC, 5hmC natively",
                       "Yes, chromosome-scale","Real-time (adaptive sampling)"),
  stringsAsFactors=FALSE
)
kable(lr_apps, caption="Long-read vs short-read capability comparison") %>%
  kable_styling(bootstrap_options=c("striped","hover","condensed"), full_width=TRUE, font_size=11) %>%
  column_spec(1, bold=TRUE)
Long-read vs short-read capability comparison
Application Short_Read PacBio_HiFi Oxford_Nanopore
Structural variant calling Limited (Manta: >50bp only) Excellent (PBSV, Sniffles2) Excellent (Sniffles2)
Full-length isoform sequencing (Iso-Seq) No (splice graphs only) Yes (Iso-Seq 3) Yes (IsoSeq-like)
Direct RNA sequencing No (cDNA proxy) No Yes (direct RNA)
De novo genome assembly No (fragmented assemblies) Best accuracy (hifiasm) Good (Flye, Verkko)
Repeat expansion genotyping Partial (ExpansionHunter) Yes (TRGT) Yes (Straglr)
CpG methylation (direct) No (bisulfite proxy only) 5mC with pb-CpG-tools 6mA, 5mC, 5hmC natively
Haplotype phasing Limited (<10kb haplotypes) Yes, chromosome-scale Yes, chromosome-scale
Metagenomics (species ID) Species-level (16S) Yes (HiFi metagenomics) Real-time (adaptive sampling)

10.2 Nanopore Basecalling and Methylation

# ---- Nanopore: basecalling with Dorado (Oxford Nanopore's current basecaller) ----
dorado basecaller \
  hac \
  pod5_pass/ \
  --reference GRCh38.fa \
  --modified-bases 5mCG_5hmCG \
  > basecalled.bam

# Extract 5mCG methylation calls
dorado summary basecalled.bam

# Modbam2bed: convert modification BAM to bedMethyl
modkit pileup \
  basecalled.bam \
  methylation.bed \
  --cpg \
  --ref GRCh38.fa \
  --threads 16

# ---- Hifiasm: PacBio HiFi de novo assembly ----
hifiasm \
  -o sample.asm \
  -t 32 \
  --hg-size 3.1g \
  sample.hifi.fastq.gz
# Output: .hap1.p_ctg.gfa, .hap2.p_ctg.gfa (haplotype-resolved)

# Convert to FASTA
awk '/^S/{print ">"$2; print $3}' sample.asm.hap1.p_ctg.gfa > hap1.fa

11 PART 9 — Multi-omics Integration

11.1 Multi-omics Data Types and Integration Strategies

multiomics <- data.frame(
  Modality = c("WGS/WES","RNA-seq","scRNA-seq","scATAC-seq","ChIP-seq",
                "ATAC-seq","Bisulfite-seq (WGBS)","Hi-C","Proteomics (TMT/LFQ)",
                "Spatial transcriptomics (Visium)","10x Multiome (RNA+ATAC)","CITE-seq"),
  Measures = c("Germline/somatic DNA variants","Bulk gene expression","Single-cell expression",
               "Chromatin accessibility per cell","Histone modifications/TF binding",
               "Open chromatin (bulk)","DNA methylation (genome-wide)","3D genome organization",
               "Protein abundance","Spatial expression + histology","RNA+ATAC in same cell",
               "Surface protein + transcriptome"),
  Key_Tool = c("GATK/DeepVariant","DESeq2/edgeR","Scanpy/Seurat","ArchR/Signac",
               "MACS3/deeptools","MACS3","Bismark/BSMAP","HiC-Pro/Juicer",
               "MaxQuant/Perseus","Squidpy/Seurat","Signac/ArchR","Totalvi/WNN"),
  Integration_Method = c("eQTL, GWAS colocalization","GSEA, pathway","Trajectory, cell-cell comm",
                           "Peak-gene linking, motif enrichment","Chromatin state segmentation",
                           "GRN inference","DMR→DEG correlation","TADs, compartments, loops",
                           "pQTL, proteogenomics","Cell-niche analysis","Joint embedding",
                           "Weighted nearest neighbour (WNN)"),
  stringsAsFactors=FALSE
)
kable(multiomics, caption="Multi-omics modalities and integration approaches") %>%
  kable_styling(bootstrap_options=c("striped","hover","condensed"), full_width=TRUE, font_size=11)
Multi-omics modalities and integration approaches
Modality Measures Key_Tool Integration_Method
WGS/WES Germline/somatic DNA variants GATK/DeepVariant eQTL, GWAS colocalization
RNA-seq Bulk gene expression DESeq2/edgeR GSEA, pathway
scRNA-seq Single-cell expression Scanpy/Seurat Trajectory, cell-cell comm
scATAC-seq Chromatin accessibility per cell ArchR/Signac Peak-gene linking, motif enrichment
ChIP-seq Histone modifications/TF binding MACS3/deeptools Chromatin state segmentation
ATAC-seq Open chromatin (bulk) MACS3 GRN inference
Bisulfite-seq (WGBS) DNA methylation (genome-wide) Bismark/BSMAP DMR→DEG correlation
Hi-C 3D genome organization HiC-Pro/Juicer TADs, compartments, loops
Proteomics (TMT/LFQ) Protein abundance MaxQuant/Perseus pQTL, proteogenomics
Spatial transcriptomics (Visium) Spatial expression + histology Squidpy/Seurat Cell-niche analysis
10x Multiome (RNA+ATAC) RNA+ATAC in same cell Signac/ArchR Joint embedding
CITE-seq Surface protein + transcriptome Totalvi/WNN Weighted nearest neighbour (WNN)

12 PART 10 — Statistics Every Bioinformatician Must Know

12.1 Multiple Testing Correction

# Visualize the multiple testing problem and BH correction
set.seed(42)
n_tests <- 20000
true_null <- 19500
true_positives <- 500

p_null <- runif(true_null)               # null p-values: uniform
p_alt  <- rbeta(true_positives, 0.3, 8) # alternative: enriched near 0

all_p <- c(p_null, p_alt)
is_true_pos <- c(rep(FALSE, true_null), rep(TRUE, true_positives))

# BH correction
p_adj_bh <- p.adjust(all_p, method="BH")
# Bonferroni correction
p_adj_bon <- p.adjust(all_p, method="bonferroni")

results <- data.frame(
  p_raw = all_p,
  p_bh = p_adj_bh,
  p_bon = p_adj_bon,
  true_positive = is_true_pos
)

# Compare thresholds
thresh <- 0.05
cat(sprintf("At raw p < 0.05:\n  True Positives called:   %d\n  False Positives called:  %d\n  FDR: %.1f%%\n",
    sum(all_p < thresh & is_true_pos), sum(all_p < thresh & !is_true_pos),
    100*sum(all_p < thresh & !is_true_pos)/sum(all_p < thresh)))
#> At raw p < 0.05:
#>   True Positives called:   379
#>   False Positives called:  1030
#>   FDR: 73.1%
cat(sprintf("\nAt BH padj < 0.05:\n  True Positives called:   %d\n  False Positives called:  %d\n",
    sum(p_adj_bh < thresh & is_true_pos), sum(p_adj_bh < thresh & !is_true_pos)))
#> 
#> At BH padj < 0.05:
#>   True Positives called:   78
#>   False Positives called:  2
cat(sprintf("\nAt Bonferroni p < 0.05:\n  True Positives called:   %d\n  False Positives called:  %d\n",
    sum(p_adj_bon < thresh & is_true_pos), sum(p_adj_bon < thresh & !is_true_pos)))
#> 
#> At Bonferroni p < 0.05:
#>   True Positives called:   22
#>   False Positives called:  0
# P-value histogram (should be uniform under null + spike at low values for true positives)
hist_df <- data.frame(p = all_p, Source = ifelse(is_true_pos, "True Positive", "True Null"))
ggplot(hist_df, aes(x=p, fill=Source)) +
  geom_histogram(bins=50, position="stack", alpha=0.8) +
  scale_fill_manual(values=c("True Positive"="#e74c3c","True Null"="#95a5a6")) +
  labs(title="P-value Distribution Under H0 and H1",
       subtitle="Uniform under null; spike near 0 for true positives — classic diagnostic plot",
       x="P-value", y="Count", fill="") +
  theme_minimal(base_size=12)

Reference: Benjamini Y., Hochberg Y. (1995). JRSS-B 57(1):289–300. [Link]

12.2 Negative Binomial vs Poisson — Why It Matters

# Demonstrate overdispersion in RNA-seq counts
set.seed(42)
n_samp <- 1000
mu <- 100  # mean expression

poi_counts <- rpois(n_samp, mu)
nb_counts  <- rnbinom(n_samp, mu=mu, size=2)  # size = 1/dispersion; size=2 = moderate overdispersion

plot_df <- data.frame(
  Count = c(poi_counts, nb_counts),
  Distribution = rep(c("Poisson (var=mean)","Negative Binomial (var>mean)"), each=n_samp)
)

ggplot(plot_df, aes(x=Count, fill=Distribution)) +
  geom_histogram(bins=60, position="identity", alpha=0.6) +
  geom_vline(xintercept=mu, linetype="dashed", color="black", linewidth=1) +
  scale_fill_manual(values=c("Poisson (var=mean)"="#3498db",
                              "Negative Binomial (var>mean)"="#e74c3c")) +
  annotate("text", x=150, y=35, 
           label=sprintf("Poisson var: %.0f\nNB var: %.0f", var(poi_counts), var(nb_counts)),
           size=3.5) +
  labs(title="Poisson vs. Negative Binomial Count Distributions (μ=100)",
       subtitle="RNA-seq counts are overdispersed → NB model required",
       x="Count", y="Frequency", fill="") +
  theme_minimal(base_size=12) + theme(legend.position="bottom")

13 PART 11 — Industry Interview Questions (Comprehensive)

13.1 Questions from Top Genomics & Bioinformatics Companies

Compiled from Glassdoor, compbiojobs.com, Levels.fyi, LinkedIn posts, and direct reports from interviewees at: Illumina, 10x Genomics, Genentech/Roche, AstraZeneca, Moderna, Novartis, Broad Institute, DNAnexus, Sema4, Blueprint Genetics, GeneDx, and academic genomics centers.

13.1.1 Foundational Statistics & Methods

Q1 [Genentech/Roche, AZ]: “You ran 20,000 differential expression tests. A colleague filters at p<0.05 and gets 1,800 genes. What’s wrong and how do you fix it?”

Model Answer: At α=0.05, you expect 20,000 × 0.05 = 1,000 false positives by random chance — the multiple testing problem. The 1,800 “significant” genes are probably ~1,000 real + ~800 from noise (not far off). Fix: apply Benjamini-Hochberg FDR correction (DESeq2 outputs padj). BH controls the expected proportion of false positives among your called significant genes at the chosen threshold. Standard is padj < 0.05 (5% FDR). For fewer false positives, use 0.01 or 0.001. Bonferroni is overly conservative for genomics and reduces power unnecessarily. Also check the p-value histogram — under a good experiment it should be uniform + spike near 0.

Q2 [Illumina Software, DNAnexus]: “Explain the difference between sensitivity and specificity, and PPV and NPV. Why does PPV matter more in rare-disease diagnostics?”

Model Answer: Sensitivity = TP/(TP+FN) — how often the test correctly identifies true positives. Specificity = TN/(TN+FP) — how often it correctly identifies true negatives. These are intrinsic to the test. PPV (Positive Predictive Value) = TP/(TP+FP) — how often a positive result is a true positive. NPV = TN/(TN+FN). PPV and NPV depend on prevalence. For a rare disease (1-in-10,000), even a highly specific test (99.9%) gives many false positives (PPV ~9%) because the population is overwhelmingly unaffected. This is Bayes’ theorem applied clinically — the reason clinical labs report variants with population frequency context and use pre-test probability.

Q3 [Blueprint Genetics, GeneDx]: “A BRCA1 variant has PM2 (absent in gnomAD), PP3 (damaging in silico), and PVS1 (frameshift in LOF-established gene). Classify it.”

Model Answer: Per ACMG/AMP rules: PVS1=Very Strong Pathogenic, PM2=Moderate Pathogenic, PP3=Supporting Pathogenic. Combining: PVS1 + 1 Moderate + 1 Supporting = Pathogenic (reaches the 1 Very Strong + ≥1 Moderate OR 1 Very Strong + ≥2 Supporting threshold). This would be reported as Pathogenic, with recommendation for genetic counselling and cascade testing of family members. Note: PVS1 can only be applied if LOF is a known disease mechanism for the gene — BRCA1 qualifies.

Q4 [Broad Institute, NYGC]: “What is VQSR and when would you use hard filtering instead?”

Model Answer: VQSR (Variant Quality Score Recalibration) trains a Gaussian mixture model on annotations from known-true variants (HapMap, 1000G, dbSNP) to learn what “true variants” look like in the multidimensional annotation space (QD, FS, MQ, ReadPosRankSum, SOR). It scores each candidate by how true-like it looks, then tranches the callset by truth sensitivity (e.g., “99.5% of HapMap sites retained”). VQSR is preferred for cohorts ≥30 samples (needs enough data to train the model robustly). Hard filtering (single-dimensional thresholds per annotation) is used for small cohorts (<30 samples), single samples, or non-human organisms where no equivalent training resource exists. Never use VQSR on whole-exome data for SNPs without at least 5,000+ samples.

13.1.2 Pipeline Design and Architecture

Q5 [10x Genomics, Sema4]: “A scRNA-seq experiment shows two distinct UMAP clusters that you suspect is a batch effect, not biology. How do you distinguish and correct it?”

Model Answer: Distinguishing: (1) color the UMAP by batch variable — if clusters are batch-pure, it’s likely technical; (2) check whether canonical cell-type markers are distributed across clusters (if T cell markers appear in both clusters, they’re likely the same cell type split by batch); (3) run kBET or a mixing score — low mixing = batch effect; (4) check if cluster-defining genes in DE analysis are housekeeping/technical (ribosomal, mitochondrial) rather than biological.
Correction: Harmony iteratively adjusts PCA embeddings to remove batch variance (best for transcript-level batch); Seurat CCA/RPCA uses anchor-based integration; scVI (variational autoencoder) models batch as a latent variable and corrects at the count level. Post-integration QC: each cluster should contain cells from all batches mixed.

Q6 [Illumina DRAGEN, Google Genomics]: “What is the algorithmic advantage of BWA-MEM’s FM-index over a naive string search? Give the time complexity.”

Model Answer: Naive exact search: O(genome_length × read_length) per read — ~3×10⁹ × 150 = 4.5×10¹¹ operations per read, infeasible at scale. BWA’s FM-index (Burrows-Wheeler Transform + compressed suffix array): backward search narrows a range in the suffix array one query character at a time, giving O(read_length) = O(150) exact match search. The BWT has an additional advantage: it compresses the genome to ~25-35% of its original size while allowing O(m) search. For approximate matching (mismatches + indels), BWA-MEM seeds with exact matches then extends with Smith-Waterman: O(m²) per seed extension, but with heuristic limits making it practically O(m) per read on typical data.

Q7 [AstraZeneca, Novartis Computational Biology]: “Design a variant-calling pipeline from scratch for a 100-patient oncology WES study. What are the key decisions and quality checkpoints?”

Model Answer: Key decisions: (1) Reference: GRCh38 (clinical standard); (2) Aligner: BWA-MEM2 (speed) or DRAGEN (accuracy); (3) Germline caller: HaplotypeCaller in GVCF + joint calling (captures all 100 patients jointly); (4) Somatic caller: Mutect2 with tumor-normal pairs + Panel of Normals; (5) SV: Manta for short reads; (6) Annotation: VEP + gnomAD + ClinVar + CADD + SpliceAI.
Quality checkpoints: FASTQ QC (FastQC/MultiQC) → alignment QC (flagstat, mosdepth, insert size, on-target fraction) → variant QC (Ti/Tv ratio, callset size, dbSNP concordance) → biological QC (PCA of samples — expect ethnicity/ancestry clustering, flag outliers) → clinical QC (ACMG classification, variant review) → concordance with prior assay if available (positive controls).
Orchestration: Nextflow/nf-core/sarek for reproducibility; containerized (Docker/Singularity); documented in git; all software versions in pipeline_info for publication compliance.

Q8 [Genentech Bioinfo, CZ Biohub]: “What is cell-cell communication inference and what are its limitations?”

Model Answer: Tools like CellChat and CellPhoneDB cross-reference the expression of curated ligand-receptor pairs (from known databases) against scRNA-seq expression profiles. For each ordered cell-type pair (sender→receiver), they test whether the ligand is significantly expressed in the sender cluster and the receptor significantly expressed in the receiver cluster — often via permutation tests shuffling cell-type labels. Statistical significance indicates the pair could be communicating via that interaction.
Critical limitations: (1) It is a correlative hypothesis, not proof of active signalling — spatial proximity, protein secretion, and functional validation are all required; (2) the ligand-receptor database is incomplete and species-transferred; (3) it ignores spatial context (cells must actually be adjacent to signal — single-cell data doesn’t tell you this); (4) single-cell transcriptomics measures mRNA abundance, not protein levels or secretion rates; (5) results are sensitive to clustering resolution and normalization.

13.1.3 Coding and Algorithmic Questions

Q9 [Broad Institute coding screen, DNAnexus]: “Write a function to parse a CIGAR string and compute the number of reference bases consumed.”

def cigar_ref_length(cigar_string):
    """
    Compute the number of reference bases consumed by a CIGAR string.
    Reference-consuming operations: M, D, N, =, X (not I, S, H, P)
    
    Example: 50M2D10M3I40M
    -> 50(M) + 2(D) + 10(M) + 0(I-skipped) + 40(M) = 102 ref bases
    
    Time complexity: O(len(CIGAR_string))
    Space complexity: O(1)
    """
    import re
    ref_ops = set('MDN=X')  # operations that consume reference
    total = 0
    for length, op in re.findall(r'(\d+)([MIDNSHP=X])', cigar_string):
        if op in ref_ops:
            total += int(length)
    return total

# Tests
tests = [
    ("76M",       76),
    ("50M2D10M",  62),
    ("10S80M5S",  80),   # soft-clips don't consume reference
    ("40M3I37M",  77),   # insertions don't consume reference
    ("100M",     100),
]
for cigar, expected in tests:
    result = cigar_ref_length(cigar)
    status = "✓" if result == expected else "✗"
    print(f"  {status} cigar_ref_length('{cigar}') = {result} (expected {expected})")

Q10 [General bioinformatics screen]: “Given a list of genomic intervals, find all overlapping pairs efficiently.”

def find_overlapping_intervals(intervals):
    """
    Find all pairs of overlapping intervals using a sweep line algorithm.
    Time: O(n log n) — dominated by sorting
    Space: O(n) for active set
    
    Input: list of (chrom, start, end, name) tuples
    Returns: list of overlapping pairs
    """
    from sortedcontainers import SortedList  # pip install sortedcontainers
    
    # Create events: (pos, type, interval_index)
    # type: 0=start, 1=end (process starts before ends at same position)
    events = []
    for i, (chrom, start, end, name) in enumerate(intervals):
        events.append((chrom, start, 0, i))   # start event
        events.append((chrom, end,   1, i))   # end event
    events.sort()
    
    active = set()
    overlaps = []
    
    for chrom, pos, etype, idx in events:
        if etype == 0:  # interval starts
            for active_idx in active:
                a_chrom, a_start, a_end, a_name = intervals[active_idx]
                if a_chrom == chrom:  # same chromosome
                    overlaps.append((intervals[idx], intervals[active_idx]))
            active.add(idx)
        else:           # interval ends
            active.discard(idx)
    
    return overlaps

# Example
intervals = [
    ("chr1", 100, 200, "gene_A"),
    ("chr1", 150, 300, "gene_B"),  # overlaps A
    ("chr1", 250, 400, "gene_C"),  # overlaps B
    ("chr2", 100, 200, "gene_D"),  # different chrom, no overlap
]

overlaps = find_overlapping_intervals(intervals)
print(f"Found {len(overlaps)} overlapping pairs:")
for a, b in overlaps:
    print(f"  {a[3]} ({a[1]}-{a[2]}) ∩ {b[3]} ({b[1]}-{b[2]})")

14 PART 12 — Workflow Orchestration: Nextflow DSL2

14.1 Complete, Buildable Pipeline

# ---- Install Nextflow ----
curl -s https://get.nextflow.io | bash
chmod +x nextflow && sudo mv nextflow /usr/local/bin/
nextflow -version  # should be ≥24.x, DSL2 default

# ---- Install Docker ----
# https://docs.docker.com/engine/install/

# ---- Run nf-core/sarek (production WGS variant calling) ----
nextflow run nf-core/sarek \
  -revision 3.4.0 \
  -profile docker \
  --input samplesheet.csv \
  --genome GATK.GRCh38 \
  --tools haplotypecaller,vep,snpeff \
  --outdir results/ \
  -resume \
  -with-report report.html \
  -with-trace trace.txt \
  -with-dag flowchart.png

14.1.1 Nextflow DSL2 Module Example

// modules/gatk_haplotypecaller.nf
process GATK_HAPLOTYPECALLER {
    tag         "${meta.id}"
    label       'process_high'
    container   'broadinstitute/gatk:4.5.0.0'
    publishDir  "${params.outdir}/gvcf/${meta.id}", mode: 'copy'

    input:
    tuple val(meta), path(bam), path(bai)
    path  fasta
    path  fai
    path  dict
    path  dbsnp
    path  dbsnp_tbi

    output:
    tuple val(meta), path("*.g.vcf.gz"), path("*.g.vcf.gz.tbi"), emit: gvcf
    path  "versions.yml",                                          emit: versions

    script:
    """
    gatk HaplotypeCaller \\
        -R ${fasta} \\
        -I ${bam} \\
        -O ${meta.id}.g.vcf.gz \\
        -ERC GVCF \\
        --dbsnp ${dbsnp} \\
        --native-pair-hmm-threads ${task.cpus} \\
        --tmp-dir .

    cat <<-END_VERSIONS > versions.yml
    "${task.process}":
        gatk4: \$(echo \$(gatk --version 2>&1) | sed 's/^.*(GATK) v//; s/ .*\$//')
    END_VERSIONS
    """
}

Reference: Di Tommaso P. et al. (2017). Nature Biotechnology 35:316–319; Ewels P.A. et al. (2020). Nature Biotechnology 38:276–278. [nf-core]

15 PART 13 — Clinical Genomics and Regulatory Framework

15.1 CLIA/CAP Requirements for Clinical Sequencing

Clinical NGS laboratories in the United States must be:
- CLIA-certified (Clinical Laboratory Improvement Amendments) — regulates analytical validity
- CAP-accredited (College of American Pathologists) — rigorous proficiency testing + inspection
- NY DOH approved (if reporting to NY patients) — separate state requirement
- Adhere to AMP/ACMG/CAP guidelines for NGS validation: analytical sensitivity ≥99%, analytical specificity ≥99.99% for reportable positions; minimum coverage ≥20× at all reportable positions (commonly ≥250× for clinical WES)

Reference: Gargis A.S. et al. (2012). “Assuring the quality of next-generation sequencing in clinical laboratory practice.” Nature Biotechnology 30:1033–1036.

15.2 Variant Reporting Tiers (Somatic — CGC/VICC)

tiers <- data.frame(
  Tier = c("Tier I","Tier II","Tier III","Tier IV"),
  Level_of_Evidence = c(
    "Strong clinical significance — FDA-approved therapy, included in clinical guidelines",
    "Potential clinical significance — investigational therapy, off-label evidence, clinical trials",
    "Unknown significance — variants in cancer genes without strong evidence",
    "Benign or likely benign somatic variation"
  ),
  Examples = c(
    "EGFR L858R → erlotinib; BRAF V600E → vemurafenib; BRCA1/2 pathogenic → PARP inhibitor",
    "ATM loss-of-function (PARP inhibitor exploratory); TP53 GOF (emerging prognostic)",
    "Novel missense in ARID1A; VUS in POLE",
    "Synonymous variants; common population polymorphisms in cancer genes"
  ),
  stringsAsFactors=FALSE
)
kable(tiers, caption="CGC/VICC somatic variant reporting tiers (Li et al. 2017)") %>%
  kable_styling(bootstrap_options=c("striped","hover"), full_width=TRUE) %>%
  column_spec(1, bold=TRUE, color="white",
              background=c("#e74c3c","#e67e22","#f39c12","#27ae60"))
CGC/VICC somatic variant reporting tiers (Li et al. 2017)
Tier Level_of_Evidence Examples
Tier I Strong clinical significance — FDA-approved therapy, included in clinical guidelines EGFR L858R → erlotinib; BRAF V600E → vemurafenib; BRCA1/2 pathogenic → PARP inhibitor
Tier II Potential clinical significance — investigational therapy, off-label evidence, clinical trials ATM loss-of-function (PARP inhibitor exploratory); TP53 GOF (emerging prognostic)
Tier III Unknown significance — variants in cancer genes without strong evidence Novel missense in ARID1A; VUS in POLE
Tier IV Benign or likely benign somatic variation Synonymous variants; common population polymorphisms in cancer genes

Reference: Li M.M. et al. (2017). “Standards and Guidelines for the Interpretation and Reporting of Sequence Variants in Cancer.” JMD 19(1):4–23.

16 PART 14 — File Format Master Reference

formats <- data.frame(
  Format = c(".fastq.gz",".sam",".bam",".cram",".g.vcf.gz",".vcf.gz",".bcf",
             ".maf",".gtf/.gff3",".bed",".mtx",".h5/.h5ad",".rds",
             ".bigwig/.bw",".bedgraph",".narrowPeak",".broadPeak"),
  Full_Name = c("FASTQ (gzip-compressed)","Sequence Alignment/Map (text)","Binary Alignment/Map",
                 "CRAM (reference-compressed)","Genomic VCF","Variant Call Format (bgzipped)",
                 "Binary VCF","Mutation Annotation Format","Gene Transfer/General Feature Format",
                 "Browser Extensible Data","MatrixMarket sparse","HDF5 / AnnData",
                 "R Data Serialization","BigWig binary track","BEDgraph signal track",
                 "Narrow peak calls","Broad peak calls"),
  Stage = c("Raw/trimmed reads","Aligned reads (text)","Aligned reads (binary)",
             "Aligned reads (space-efficient)","Per-sample variant evidence",
             "Variant calls (indexed)","Binary VCF (indexed)","Cancer cohort variants",
             "Gene annotation","Genomic intervals","scRNA-seq counts (sparse)",
             "Analysis objects (scRNA-seq)","Seurat/R objects","Normalized signal (ATAC, ChIP)",
             "Normalized signal (text)","ChIP/ATAC peak calls","ChIP broad domains"),
  Key_Tool = c("fastp,FastQC","bwa-mem2,minimap2","samtools","samtools view -C",
               "GATK HC","GATK,FreeBayes","bcftools","vcf2maf","GENCODE,Ensembl",
               "BEDtools","STARsolo,Cell Ranger","Scanpy,AnnData","Seurat",
               "deeptools bamCoverage","deeptools","MACS3","MACS3 --broad"),
  Index = c("none (gz stream)","none",".bai or .csi",".crai",".tbi",".tbi",".csi",
            "none",".gtf.gz + .tbi (optional)",".bai (AILIST)",
            "features.tsv + barcodes.tsv","self-indexed","none",".bw self-indexed","none","none","none"),
  stringsAsFactors=FALSE
)
kable(formats, caption="Complete genomics file format reference") %>%
  kable_styling(bootstrap_options=c("striped","hover","condensed"), full_width=TRUE, font_size=11) %>%
  column_spec(1, bold=TRUE, color="white", background="#2c3e50", monospace=TRUE)
Complete genomics file format reference
Format Full_Name Stage Key_Tool Index
.fastq.gz FASTQ (gzip-compressed) Raw/trimmed reads fastp,FastQC none (gz stream)
.sam Sequence Alignment/Map (text) Aligned reads (text) bwa-mem2,minimap2 none
.bam Binary Alignment/Map Aligned reads (binary) samtools .bai or .csi
.cram CRAM (reference-compressed) Aligned reads (space-efficient) samtools view -C .crai
.g.vcf.gz Genomic VCF Per-sample variant evidence GATK HC .tbi
.vcf.gz Variant Call Format (bgzipped) Variant calls (indexed) GATK,FreeBayes .tbi
.bcf Binary VCF Binary VCF (indexed) bcftools .csi
.maf Mutation Annotation Format Cancer cohort variants vcf2maf none
.gtf/.gff3 Gene Transfer/General Feature Format Gene annotation GENCODE,Ensembl .gtf.gz + .tbi (optional)
.bed Browser Extensible Data Genomic intervals BEDtools .bai (AILIST)
.mtx MatrixMarket sparse scRNA-seq counts (sparse) STARsolo,Cell Ranger features.tsv + barcodes.tsv
.h5/.h5ad HDF5 / AnnData Analysis objects (scRNA-seq) Scanpy,AnnData self-indexed
.rds R Data Serialization Seurat/R objects Seurat none
.bigwig/.bw BigWig binary track Normalized signal (ATAC, ChIP) deeptools bamCoverage .bw self-indexed
.bedgraph BEDgraph signal track Normalized signal (text) deeptools none
.narrowPeak Narrow peak calls ChIP/ATAC peak calls MACS3 none
.broadPeak Broad peak calls ChIP broad domains MACS3 –broad none

17 PART 15 — Key Tools Reference Table

tools <- data.frame(
  Tool = c("FastQC","MultiQC","fastp","Trim Galore","BWA-MEM2","STAR","HISAT2","minimap2",
           "samtools","GATK4","DeepVariant","Mutect2","Strelka2","Manta","Sniffles2",
           "CNVkit","VEP","SnpEff","ANNOVAR","DESeq2","edgeR","Salmon","kallisto",
           "featureCounts","Cell Ranger","STARsolo","Scanpy","Seurat","Harmony",
           "Nextflow","nf-core","Docker/Singularity","mosdepth","picard","bcftools"),
  Category = c("QC","QC","Trimming","Trimming","Alignment","Alignment","Alignment","Alignment (LR)",
               "BAM handling","Variant calling","Variant calling (DL)","Somatic calling",
               "Germline/somatic","SV calling","SV (long read)","CNV","Annotation","Annotation",
               "Annotation","DE analysis","DE analysis","Quantification","Quantification",
               "Quantification","scRNA-seq","scRNA-seq","scRNA-seq","scRNA-seq","Integration",
               "Workflow","Workflow library","Containers","Coverage","Library prep QC","VCF handling"),
  Language = c("Java","Python","C++","Python/shell","C","C","C","C","C","Java",
                "Python/C++","Java","C++","C++","C++","Python","Perl/Python","Java",
                "Perl","R/Bioc","R/Bioc","C++","C++","C","Python","C","Python","R","R",
                "Groovy","Nextflow/Python","—","C","Java","C"),
  Key_Reference = c("Andrews 2010","Ewels 2016","Chen 2018","Krueger","Vasimuddin 2019",
                    "Dobin 2013","Kim 2015","Li 2018","Li 2009","McKenna 2010",
                    "Poplin 2018 NBT","Benjamin 2019","Kim 2018","Chen 2016","Smolka 2022",
                    "Talevich 2016","McLaren 2016","Cingolani 2012","Wang 2010","Love 2014",
                    "Robinson 2010","Patro 2017","Bray 2016","Liao 2014",
                    "10x Genomics","Kaminow 2021","Wolf 2018","Hao 2021 Cell","Korsunsky 2019",
                    "Di Tommaso 2017","Ewels 2020","—","Pedersen 2018","Picard Broad","Li 2011"),
  stringsAsFactors=FALSE
)
kable(tools, caption="Comprehensive bioinformatics tool reference") %>%
  kable_styling(bootstrap_options=c("striped","hover","condensed"), full_width=TRUE, font_size=11) %>%
  column_spec(1, bold=TRUE, monospace=TRUE) %>%
  column_spec(2, color="white", 
              background=c(rep("#2980b9",2), rep("#16a085",2), rep("#8e44ad",4),
                           rep("#34495e",1), rep("#c0392b",7), rep("#e67e22",3),
                           rep("#2ecc71",2), rep("#1abc9c",4), rep("#9b59b6",4),
                           rep("#f39c12",1), rep("#e74c3c",2), rep("#95a5a6",3)))
Comprehensive bioinformatics tool reference
Tool Category Language Key_Reference
FastQC QC Java Andrews 2010
MultiQC QC Python Ewels 2016
fastp Trimming C++ Chen 2018
Trim Galore Trimming Python/shell Krueger
BWA-MEM2 Alignment C Vasimuddin 2019
STAR Alignment C Dobin 2013
HISAT2 Alignment C Kim 2015
minimap2 Alignment (LR) C Li 2018
samtools BAM handling C Li 2009
GATK4 Variant calling Java McKenna 2010
DeepVariant Variant calling (DL) Python/C++ Poplin 2018 NBT
Mutect2 Somatic calling Java Benjamin 2019
Strelka2 Germline/somatic C++ Kim 2018
Manta SV calling C++ Chen 2016
Sniffles2 SV (long read) C++ Smolka 2022
CNVkit CNV Python Talevich 2016
VEP Annotation Perl/Python McLaren 2016
SnpEff Annotation Java Cingolani 2012
ANNOVAR Annotation Perl Wang 2010
DESeq2 DE analysis R/Bioc Love 2014
edgeR DE analysis R/Bioc Robinson 2010
Salmon Quantification C++ Patro 2017
kallisto Quantification C++ Bray 2016
featureCounts Quantification C Liao 2014
Cell Ranger scRNA-seq Python 10x Genomics
STARsolo scRNA-seq C Kaminow 2021
Scanpy scRNA-seq Python Wolf 2018
Seurat scRNA-seq R Hao 2021 Cell
Harmony Integration R Korsunsky 2019
Nextflow Workflow Groovy Di Tommaso 2017
nf-core Workflow library Nextflow/Python Ewels 2020
Docker/Singularity Containers — —
mosdepth Coverage C Pedersen 2018
picard Library prep QC Java Picard Broad
bcftools VCF handling C Li 2011

18 PART 16 — THE FULL CURRICULUM: Beginner → Pre-Doctoral Computational Genomics

18.1 A Standalone Deep-Dive for Computer Science & Bioinformatics Master’s Graduates

This part is written as a self-contained second course layered underneath everything above. It starts from first principles (what is a base pair, what is Big-O, what is a Markov chain) and climbs to the level expected of an incoming PhD student or senior computational biologist — including the CS theory (data structures, complexity, probabilistic models, machine learning) that underlies every tool used in Parts 1–15. Every subsection is tiered: Beginner → Intermediate → Advanced/Pre-doctoral, each with runnable R and Python code.

18.2 16.0 How to Use This Part

library(ggplot2)

nodes <- data.frame(
  concept = c("Molecular Biology\nBasics","Probability &\nStatistics","Data Structures\n& Algorithms",
              "String Algorithms\n(BWT/Suffix Arrays)","Sequence Alignment\nTheory","Hidden Markov\nModels",
              "Graph Theory\n(De Bruijn/OLC)","Linear Algebra\n(PCA/SVD)","Machine Learning\nFoundations",
              "Deep Learning\n(CNN/Transformers)","Bayesian Statistics","Population Genetics",
              "Systems & Distributed\nComputing","Software Engineering\n& Reproducibility"),
  x = c(1,1,1, 2,2,2.5, 3,3,3.5, 4.5,4,4, 5.5,5.5),
  y = c(5,4,3, 4.5,3.5,2.5, 3.5,2,1.5, 1.2,2.8,4.2, 2.2,3.8),
  tier = c("Foundation","Foundation","Foundation","Core CS","Core CS","Core CS",
           "Core CS","Core CS","Advanced","Advanced","Advanced","Advanced","Systems","Systems")
)

edges <- data.frame(
  x=c(1,1,1,1,2,2,2.5,3,3,3.5,1,4,3,5.5),
  y=c(5,4,3,3,4.5,3.5,2.5,3.5,2,1.5,4,2.8,2,2.2),
  xend=c(2,2,2,2.5,2.5,3,3.5,3.5,3.5,4.5,4,4.5,5.5,4.5),
  yend=c(4.5,3.5,3.5,2.5,3.5,3.5,1.5,2,1.5,1.2,2.8,1.2,2.2,1.2)
)

ggplot() +
  geom_segment(data=edges, aes(x=x,y=y,xend=xend,yend=yend), color="grey70", linewidth=0.6,
               arrow=arrow(length=unit(0.15,"cm"))) +
  geom_point(data=nodes, aes(x=x,y=y,color=tier), size=14, alpha=0.85) +
  geom_text(data=nodes, aes(x=x,y=y,label=concept), size=2.6, fontface="bold", color="white") +
  scale_color_manual(values=c("Foundation"="#3498db","Core CS"="#e67e22",
                               "Advanced"="#9b59b6","Systems"="#16a085")) +
  labs(title="Concept Dependency Map: Beginner → Pre-Doctoral",
       subtitle="Arrows show 'is a prerequisite for' — follow left to right",
       color="Tier") +
  theme_void(base_size=12) + theme(legend.position="bottom", plot.title=element_text(face="bold"))
Concept dependency map for this curriculum

Concept dependency map for this curriculum

18.3 16.1 TIER 1 — Beginner Foundations

18.3.1 16.1.1 Molecular Biology, From Zero

What is a genome, precisely? A genome is a sequence over the 4-letter alphabet Σ = {A, C, G, T} (DNA) representing the complete hereditary information of an organism. The human genome is ~3.1 × 10⁹ base pairs, organized into 23 chromosome pairs (22 autosomes + X/Y). DNA is double-stranded and complementary: A pairs with T, C pairs with G (Watson-Crick base pairing), and the two strands run in opposite (antiparallel) 5′→3′ directions — this is why every aligner must consider both the forward and reverse-complement strand of a read.

Central Dogma: DNA → (transcription) → RNA → (translation) → Protein. This one-directional flow (with reverse transcription as the well-known retroviral exception) is why RNA-seq (Part 5–6) requires fundamentally different alignment strategies than DNA-seq (Part 2–3): mRNA has introns spliced out, so a read spanning an exon-exon junction has no contiguous match in genomic DNA.

Reference: Crick F. (1970). “Central dogma of molecular biology.” Nature 227:561–563. [Link]

library(ggplot2)
stages <- data.frame(
  x = c(1,2,3), y = c(1,1,1),
  label = c("DNA\n(genome)","RNA\n(transcriptome)","Protein\n(proteome)"),
  color = c("#3498db","#e74c3c","#27ae60")
)
ggplot() +
  geom_segment(data=data.frame(x=c(1.3,2.3), xend=c(1.7,2.7), y=c(1,1), yend=c(1,1)),
               aes(x=x,xend=xend,y=y,yend=yend), arrow=arrow(length=unit(0.3,"cm")), linewidth=1.2) +
  annotate("text", x=1.5, y=1.15, label="Transcription", size=3.5, fontface="italic") +
  annotate("text", x=2.5, y=1.15, label="Translation", size=3.5, fontface="italic") +
  geom_point(data=stages, aes(x=x,y=y), size=30, color=stages$color, alpha=0.85) +
  geom_text(data=stages, aes(x=x,y=y,label=label), color="white", fontface="bold", size=4) +
  xlim(0.5,3.5) + ylim(0.7,1.3) +
  labs(title="The Central Dogma of Molecular Biology") +
  theme_void(base_size=13)

Beginner exercise — reverse complement (every bioinformatician’s “hello world”):

def reverse_complement(seq: str) -> str:
    """
    Compute the reverse complement of a DNA sequence.
    A<->T, C<->G, then reverse the order.
    O(n) time, O(n) space.
    """
    complement = {'A':'T','T':'A','C':'G','G':'C','N':'N',
                  'a':'t','t':'a','c':'g','g':'c','n':'n'}
    return ''.join(complement[base] for base in reversed(seq))

seq = "ATCGGGCATNNTAC"
print(f"Original:            {seq}")
print(f"Reverse complement:  {reverse_complement(seq)}")
# ATCGGGCATNNTAC -> GTANNATGCCCGAT
reverse_complement <- function(seq) {
  comp_map <- c(A="T", T="A", C="G", G="C", N="N")
  chars <- strsplit(seq, "")[[1]]
  comp <- comp_map[chars]
  paste(rev(comp), collapse = "")
}
cat("Reverse complement of ATCGGGCAT:", reverse_complement("ATCGGGCAT"), "\n")
#> Reverse complement of ATCGGGCAT: ATGCCCGAT
# Using Biostrings (the Bioconductor standard, used throughout industry)
if (requireNamespace("Biostrings", quietly=TRUE)) {
  library(Biostrings)
  dna <- DNAString("ATCGGGCATNNTAC")
  cat("Biostrings reverseComplement():", as.character(reverseComplement(dna)), "\n")
  cat("GC content:", letterFrequency(dna, "GC", as.prob=TRUE), "\n")
}
#> Biostrings reverseComplement(): GTANNATGCCCGAT 
#> GC content: 0.4285714

18.3.2 16.1.2 Codons, Reading Frames, and the Genetic Code

codon_table <- data.frame(
  Codon = c("ATG","TAA","TAG","TGA","TTT","GGG"),
  Amino_Acid = c("Methionine (Met, M) — START codon","STOP","STOP","STOP",
                  "Phenylalanine (Phe, F)","Glycine (Gly, G)"),
  Note = c("Translation initiation site","","","",
           "Codon degeneracy: 61 sense codons encode only 20 amino acids",
           "")
)
kable(codon_table, caption="Genetic code examples — the redundancy (degeneracy) of the code matters for synonymous variant interpretation") %>%
  kable_styling(bootstrap_options="striped", full_width=FALSE)
Genetic code examples — the redundancy (degeneracy) of the code matters for synonymous variant interpretation
Codon Amino_Acid Note
ATG Methionine (Met, M) — START codon Translation initiation site
TAA STOP
TAG STOP
TGA STOP
TTT Phenylalanine (Phe, F) Codon degeneracy: 61 sense codons encode only 20 amino acids
GGG Glycine (Gly, G) 
# Beginner: translate DNA to protein (open reading frame)
genetic_code = {
    'TTT':'F','TTC':'F','TTA':'L','TTG':'L','CTT':'L','CTC':'L','CTA':'L','CTG':'L',
    'ATT':'I','ATC':'I','ATA':'I','ATG':'M','GTT':'V','GTC':'V','GTA':'V','GTG':'V',
    'TCT':'S','TCC':'S','TCA':'S','TCG':'S','CCT':'P','CCC':'P','CCA':'P','CCG':'P',
    'ACT':'T','ACC':'T','ACA':'T','ACG':'T','GCT':'A','GCC':'A','GCA':'A','GCG':'A',
    'TAT':'Y','TAC':'Y','TAA':'*','TAG':'*','CAT':'H','CAC':'H','CAA':'Q','CAG':'Q',
    'AAT':'N','AAC':'N','AAA':'K','AAG':'K','GAT':'D','GAC':'D','GAA':'E','GAG':'E',
    'TGT':'C','TGC':'C','TGA':'*','TGG':'W','CGT':'R','CGC':'R','CGA':'R','CGG':'R',
    'AGT':'S','AGC':'S','AGA':'R','AGG':'R','GGT':'G','GGC':'G','GGA':'G','GGG':'G',
}

def translate(seq: str) -> str:
    """Translate DNA sequence to protein, starting at position 0. O(n)."""
    protein = []
    for i in range(0, len(seq) - 2, 3):
        codon = seq[i:i+3]
        aa = genetic_code.get(codon, 'X')
        if aa == '*':
            break
        protein.append(aa)
    return ''.join(protein)

dna = "ATGGCTGATTTTGGGTAA"
print(f"DNA:     {dna}")
print(f"Protein: {translate(dna)}")   # MADFG

# Why frameshift indels are so damaging (PVS1 in ACMG!):
frameshift_dna = "ATGGCTGA" + "T" + "TTGGGTAA"  # 1bp insertion shifts everything downstream
print(f"\nWith 1bp insertion (frameshift): {translate(frameshift_dna)}")
print("-> completely different protein from the insertion point onward")

18.3.3 16.1.3 Big-O Notation for Biologists (CS Foundation)

Every algorithm in this document has a time complexity and space complexity. For a CS/bioinformatics graduate, being able to state and justify these is a baseline interview expectation.

Notation Name Genomics Example
O(1) Constant Hash table lookup of a k-mer
O(log n) Logarithmic Binary search in a sorted BAM index
O(n) Linear Streaming through a FASTQ file once
O(n log n) Linearithmic Sorting reads by coordinate (samtools sort)
O(n·m) Bilinear Naive alignment of read (m) against genome (n)
O(n²) Quadratic Naive dynamic programming alignment matrix
O(2ⁿ) Exponential Naive haplotype phasing over n heterozygous sites (brute force) 
import time
import random

def naive_search(genome: str, pattern: str) -> list:
    """O(n*m) naive substring search — the baseline every aligner improves upon."""
    positions = []
    n, m = len(genome), len(pattern)
    for i in range(n - m + 1):
        if genome[i:i+m] == pattern:
            positions.append(i)
    return positions

def benchmark_naive_search():
    random.seed(42)
    genome = ''.join(random.choices('ACGT', k=100000))
    pattern = genome[50000:50020]  # 20bp pattern known to exist
    
    start = time.time()
    hits = naive_search(genome, pattern)
    elapsed = time.time() - start
    print(f"Naive O(n*m) search: {elapsed*1000:.2f}ms for {len(genome):,}bp genome")
    print(f"Extrapolated to human genome (3.1Gb): ~{elapsed * (3.1e9/1e5):.1f}s PER READ")
    print(f"For 1 billion reads: ~{elapsed * (3.1e9/1e5) * 1e9 / 3600 / 24 / 365:.0f} YEARS")
    print("--> This is exactly why BWA's FM-index (O(m) per read) was necessary.")

benchmark_naive_search()

18.4 16.2 TIER 2 — Core Computer Science for Genomics

18.4.1 16.2.1 Suffix Arrays, Suffix Trees, and the FM-Index — Full Construction

Pre-doctoral depth: Part 2 introduced the use of the BWT/FM-index. Here we build one from scratch and analyze its complexity rigorously — this is standard technical-interview and qualifying-exam material at genomics-focused CS programs.

A suffix array SA[0..n-1] of string T (length n) is a permutation of {0,…,n-1} such that T[SA[i]:] is the i-th smallest suffix lexicographically. Naive construction: O(n² log n) (sort n suffixes, each comparison O(n)). Prefix doubling (Manber-Myers, 1993) achieves O(n log n). DC3/skew algorithm (Kärkkäinen-Sanders, 2003) achieves O(n) — linear time, used in production genome indexers.

def build_suffix_array_naive(s: str) -> list:
    """
    Naive suffix array construction: O(n^2 log n)
    Educational only — production tools use SA-IS or DC3 (O(n))
    """
    s = s + '$'  # sentinel, lexicographically smallest
    n = len(s)
    suffixes = sorted(range(n), key=lambda i: s[i:])
    return suffixes

def build_bwt_from_sa(s: str, sa: list) -> str:
    """BWT[i] = s[SA[i]-1] (the character preceding each sorted suffix)."""
    s = s + '$'
    return ''.join(s[i-1] for i in sa)

def build_fm_index(bwt: str):
    """
    Build FM-index components: C (first-occurrence array) and Occ (rank array).
    Enables O(1) rank queries -> O(m) backward search total.
    """
    chars = sorted(set(bwt))
    counts = {c: bwt.count(c) for c in chars}
    
    # C[c] = number of characters in BWT lexicographically smaller than c
    C = {}
    total = 0
    for c in sorted(chars):
        C[c] = total
        total += counts[c]
    
    # Occ[c][i] = number of occurrences of c in bwt[0:i]
    Occ = {c: [0]*(len(bwt)+1) for c in chars}
    for i, ch in enumerate(bwt):
        for c in chars:
            Occ[c][i+1] = Occ[c][i] + (1 if ch == c else 0)
    
    return C, Occ

def fm_backward_search(bwt: str, C: dict, Occ: dict, pattern: str):
    """
    Exact pattern search via FM-index backward search.
    Time: O(m) where m = len(pattern) -- independent of genome length n!
    """
    top, bot = 0, len(bwt)
    for c in reversed(pattern):
        if c not in C:
            return 0, 0
        top = C[c] + Occ[c][top]
        bot = C[c] + Occ[c][bot]
        if top >= bot:
            return 0, 0
    return top, bot

# ---- Demo ----
genome = "GATTACAGATTACAGATTACA"
sa = build_suffix_array_naive(genome)
bwt = build_bwt_from_sa(genome, sa)
C, Occ = build_fm_index(bwt)

print(f"Genome:  {genome}")
print(f"BWT:     {bwt}")
print(f"Suffix Array: {sa}\n")

for pattern in ["GATTACA", "TAC", "XYZ"]:
    top, bot = fm_backward_search(bwt, C, Occ, pattern)
    n_hits = bot - top
    print(f"Pattern '{pattern}': {n_hits} occurrence(s) [SA range {top}:{bot}]")

Reference: Ferragina P., Manzini G. (2000). “Opportunistic Data Structures with Applications.” FOCS 2000 (defines the FM-index formally). Kärkkäinen J., Sanders P. (2003). “Simple Linear Work Suffix Array Construction.” ICALP 2003. [FM-index paper]

18.4.2 16.2.2 Dynamic Programming — Needleman-Wunsch and Smith-Waterman From Scratch

Every “seed-and-extend” aligner (BWA-MEM, minimap2) uses a local alignment DP step after seeding. Understanding the recurrence is essential.

Needleman-Wunsch (global alignment, 1970): \[F(i,j) = \max \begin{cases} F(i-1,j-1) + s(x_i,y_j) & \text{(match/mismatch)} \\ F(i-1,j) - d & \text{(deletion)} \\ F(i,j-1) - d & \text{(insertion)} \end{cases}\]

Smith-Waterman (local alignment, 1981) adds a floor of 0 — alignment can restart anywhere, so a local match “wins” over extending a bad global one: \[F(i,j) = \max(0, \ldots)\]

Both are O(nm) time and O(nm) space (reducible to O(min(n,m)) space with Hirschberg’s algorithm for the traceback-free case).

import numpy as np

def smith_waterman(seq1: str, seq2: str, match=2, mismatch=-1, gap=-2):
    """
    Local alignment via Smith-Waterman dynamic programming.
    Time: O(n*m), Space: O(n*m)
    Returns: (score, aligned_seq1, aligned_seq2)
    """
    n, m = len(seq1), len(seq2)
    H = np.zeros((n+1, m+1), dtype=int)
    traceback = np.zeros((n+1, m+1), dtype=int)  # 0=stop, 1=diag, 2=up, 3=left
    
    max_score, max_pos = 0, (0, 0)
    
    for i in range(1, n+1):
        for j in range(1, m+1):
            s = match if seq1[i-1] == seq2[j-1] else mismatch
            diag = H[i-1][j-1] + s
            up   = H[i-1][j] + gap
            left = H[i][j-1] + gap
            
            H[i][j] = max(0, diag, up, left)
            
            if H[i][j] == 0:      traceback[i][j] = 0
            elif H[i][j] == diag: traceback[i][j] = 1
            elif H[i][j] == up:   traceback[i][j] = 2
            else:                 traceback[i][j] = 3
            
            if H[i][j] > max_score:
                max_score, max_pos = H[i][j], (i, j)
    
    # Traceback
    aligned1, aligned2 = [], []
    i, j = max_pos
    while traceback[i][j] != 0:
        if traceback[i][j] == 1:
            aligned1.append(seq1[i-1]); aligned2.append(seq2[j-1]); i -= 1; j -= 1
        elif traceback[i][j] == 2:
            aligned1.append(seq1[i-1]); aligned2.append('-'); i -= 1
        else:
            aligned1.append('-'); aligned2.append(seq2[j-1]); j -= 1
    
    return max_score, ''.join(reversed(aligned1)), ''.join(reversed(aligned2))

# Example: align a read against a reference with a mismatch and small indel
reference = "ACGTTGACAGGTACCA"
read      = "TTGACTGGTACC"   # note: 1 mismatch + models a small indel scenario

score, aln_ref, aln_read = smith_waterman(reference, read)
print(f"Score: {score}")
print(f"Ref:  {aln_ref}")
print(f"Read: {aln_read}")
print("(This O(n*m) DP is exactly the 'extend' step after BWA-MEM's O(m) seed step)")
# Visualize a Smith-Waterman scoring matrix as a heatmap
seq1 <- "ACGTTGAC"
seq2 <- "ACTTGAC"
n <- nchar(seq1); m <- nchar(seq2)
match_score <- 2; mismatch_score <- -1; gap_penalty <- -2

H <- matrix(0, n+1, m+1)
s1 <- strsplit(seq1, "")[[1]]
s2 <- strsplit(seq2, "")[[1]]

for (i in 2:(n+1)) {
  for (j in 2:(m+1)) {
    s <- ifelse(s1[i-1] == s2[j-1], match_score, mismatch_score)
    diag <- H[i-1, j-1] + s
    up   <- H[i-1, j] + gap_penalty
    left <- H[i, j-1] + gap_penalty
    H[i, j] <- max(0, diag, up, left)
  }
}

H_df <- as.data.frame(as.table(H))
colnames(H_df) <- c("i","j","score")
H_df$i <- as.numeric(H_df$i) - 1
H_df$j <- as.numeric(H_df$j) - 1

ggplot(H_df, aes(x=j, y=i, fill=score)) +
  geom_tile(color="white") +
  geom_text(aes(label=score), size=3) +
  scale_fill_gradient2(low="white", mid="#aed6f1", high="#e74c3c", midpoint=max(H)/2) +
  scale_y_reverse() +
  labs(title="Smith-Waterman DP Matrix", subtitle=paste0("Ref: ",seq1," | Read: ",seq2),
       x="Read position (j)", y="Reference position (i)") +
  theme_minimal(base_size=11) + theme(legend.position="none")

Reference: Needleman S.B., Wunsch C.D. (1970). J Mol Biol 48(3):443–453. Smith T.F., Waterman M.S. (1981). J Mol Biol 147(1):195–197. [Smith-Waterman]

18.4.3 16.2.3 Hidden Markov Models — the Mathematics Behind HaplotypeCaller and Phred

Pre-doctoral level: GATK’s PairHMM (used inside HaplotypeCaller, Part 3) computes the likelihood of a read given a candidate haplotype by modeling three hidden states per position: Match (M), Insertion (I), Deletion (D) — exactly the classic profile-HMM structure from Durbin et al.’s Biological Sequence Analysis (1998).

HMM formalism:
- States: \(Q = \{q_1, ..., q_N\}\) (here: M, I, D)
- Transition matrix: \(A = \{a_{ij}\} = P(q_j \text{ at } t+1 \mid q_i \text{ at } t)\)
- Emission probabilities: \(B = \{b_j(o_t)\} = P(\text{observation } o_t \mid \text{state } q_j)\)
- The Forward algorithm computes \(P(O|\lambda)\) — probability of the observed read given the model — in O(N²T) via dynamic programming (avoiding the naive O(N^T) enumeration of all state paths)

import numpy as np

def forward_algorithm(observations, states, start_prob, trans_prob, emit_prob):
    """
    HMM Forward algorithm: compute P(O|model) via dynamic programming.
    Time: O(N^2 * T) where N=#states, T=#observations
    (vs naive O(N^T) brute-force over all state sequences)
    
    This is the core computational primitive inside GATK's PairHMM
    genotype-likelihood calculation.
    """
    T = len(observations)
    N = len(states)
    alpha = np.zeros((T, N))
    
    # Initialization
    for s in range(N):
        alpha[0, s] = start_prob[s] * emit_prob[s][observations[0]]
    
    # Recursion
    for t in range(1, T):
        for s in range(N):
            alpha[t, s] = sum(alpha[t-1, s2] * trans_prob[s2][s] for s2 in range(N)) \
                          * emit_prob[s][observations[t]]
    
    # Termination
    prob = sum(alpha[T-1, s] for s in range(N))
    return prob, alpha

# Simplified 2-state HMM: "Match" vs "Mismatch/Error" emitting quality-consistent bases
states = ['Match', 'Error']
start_prob = [0.99, 0.01]
trans_prob = {0: [0.95, 0.05], 1: [0.60, 0.40]}  # errors tend to cluster (context-dependence)
# Emission: probability of observing "consistent" (0) vs "inconsistent" (1) base given state
emit_prob = {0: [0.98, 0.02], 1: [0.30, 0.70]}

observations = [0, 0, 0, 1, 0, 0]  # 0=base matches reference, 1=mismatch
prob, alpha = forward_algorithm(observations, states, start_prob, trans_prob, emit_prob)
print(f"P(observed read pattern | HMM) = {prob:.6f}")
print("\nForward table (alpha):")
print(f"{'t':>3} {'P(Match)':>12} {'P(Error)':>12}")
for t in range(len(observations)):
    print(f"{t:>3} {alpha[t,0]:>12.6f} {alpha[t,1]:>12.6f}")

Reference: Durbin R., Eddy S.R., Krogh A., Mitchison G. (1998). Biological Sequence Analysis. Cambridge University Press (the standard graduate HMM-in-bioinformatics textbook). Rabiner L.R. (1989). “A tutorial on hidden Markov models.” Proceedings of the IEEE 77(2):257–286. [Rabiner tutorial]

18.4.4 16.2.4 De Bruijn Graphs — Assembly Theory

Why HaplotypeCaller and genome assemblers both use De Bruijn graphs: a De Bruijn graph of order k represents every k-mer in the reads as a node, with edges connecting k-mers that overlap by k-1 bases. Reconstructing the original sequence becomes an Eulerian path problem (visit every edge exactly once) — solvable in O(E) time (Hierholzer’s algorithm), vastly cheaper than the NP-hard Hamiltonian path problem (visit every node once) that overlap-layout-consensus (OLC) assemblers like Celera/hifiasm-for-lower-coverage must approximate heuristically.

from collections import defaultdict

def build_de_bruijn_graph(reads, k):
    """
    Build a De Bruijn graph from k-mers.
    Nodes = (k-1)-mers; Edges = k-mers (connecting the prefix and suffix (k-1)-mers)
    Time: O(total_bases), Space: O(unique k-mers)
    """
    graph = defaultdict(list)
    for read in reads:
        for i in range(len(read) - k + 1):
            kmer = read[i:i+k]
            prefix, suffix = kmer[:-1], kmer[1:]
            graph[prefix].append(suffix)
    return graph

def find_eulerian_path(graph):
    """
    Hierholzer's algorithm for Eulerian path: O(E) time.
    This is what reconstructs the assembled sequence from the graph.
    """
    graph = {k: list(v) for k, v in graph.items()}
    
    # Find start node (out-degree - in-degree = 1, or any node with edges)
    out_deg = {n: len(v) for n, v in graph.items()}
    in_deg = defaultdict(int)
    for n, targets in graph.items():
        for t in targets:
            in_deg[t] += 1
    
    start = next(iter(graph))
    for n in graph:
        if out_deg.get(n,0) - in_deg.get(n,0) == 1:
            start = n
            break
    
    stack, path = [start], []
    while stack:
        node = stack[-1]
        if graph.get(node):
            next_node = graph[node].pop()
            stack.append(next_node)
        else:
            path.append(stack.pop())
    
    return path[::-1]

def reconstruct_sequence(path):
    """Reconstruct sequence from Eulerian path of (k-1)-mers."""
    if not path: return ""
    seq = path[0]
    for node in path[1:]:
        seq += node[-1]
    return seq

# Simulate reads from a known sequence (with overlap, as in real sequencing)
original = "ACGTACGGTCAGT"
k = 4
reads = [original[i:i+6] for i in range(0, len(original)-5, 2)]  # overlapping reads
print(f"Original sequence: {original}")
print(f"Simulated reads:   {reads}\n")

graph = build_de_bruijn_graph(reads, k)
print(f"De Bruijn graph ({k}-mers, {len(graph)} distinct {k-1}-mer nodes):")
for node, edges in list(graph.items())[:5]:
    print(f"  {node} -> {edges}")

path = find_eulerian_path(graph)
reconstructed = reconstruct_sequence(path)
print(f"\nReconstructed: {reconstructed}")
print(f"Matches original substring: {original in reconstructed or reconstructed in original}")

Reference: Pevzner P.A., Tang H., Waterman M.S. (2001). “An Eulerian path approach to DNA fragment assembly.” PNAS 98(17):9748–9753 (foundational De Bruijn graph assembly paper). [Link]
Compeau P.E.C., Pevzner P.A., Tesler G. (2011). “How to apply de Bruijn graphs to genome assembly.” Nature Biotechnology 29:987–991. [Link]

18.5 16.3 TIER 3 — Advanced Statistics and Machine Learning for Genomics

18.5.1 16.3.1 Bayesian Inference — the Foundation of Every Genotype Call

Every genotype likelihood in a VCF (Part 3, PL field) is computed via Bayes’ theorem: \[P(G \mid D) = \frac{P(D \mid G) \, P(G)}{P(D)}\] where G = genotype (0/0, 0/1, 1/1), D = observed read data. \(P(D|G)\) is the genotype likelihood (computed from base qualities via the PairHMM, §16.2.3); \(P(G)\) is the prior (typically from population allele frequency, e.g. Hardy-Weinberg expectation); \(P(D)\) is a normalizing constant. GATK reports \(-10\log_{10}P(G|D)\) as PL (Phred-scaled Likelihood) for each possible genotype, with 0 assigned to the most likely genotype.

import numpy as np
from scipy.stats import binom

def genotype_likelihood(ref_reads: int, alt_reads: int, error_rate: float = 0.01):
    """
    Compute P(D|G) for G in {0/0, 0/1, 1/1} using a binomial error model.
    This mirrors (a simplification of) what GATK's PairHMM computes.
    """
    total = ref_reads + alt_reads
    genotypes = {
        '0/0': error_rate,       # expect ~0% alt reads (all "alt" calls are errors)
        '0/1': 0.5,              # expect ~50% alt reads
        '1/1': 1 - error_rate    # expect ~100% alt reads
    }
    
    likelihoods = {}
    for gt, p_alt in genotypes.items():
        # P(observing alt_reads alt-supporting reads out of `total`, given true alt fraction p_alt)
        likelihoods[gt] = binom.pmf(alt_reads, total, p_alt)
    
    # Normalize to posterior (assuming flat prior for simplicity)
    total_lik = sum(likelihoods.values())
    posteriors = {gt: lik/total_lik for gt, lik in likelihoods.items()}
    
    # Phred-scaled likelihoods (PL), normalized so best genotype = 0
    max_lik = max(likelihoods.values())
    pl = {gt: round(-10 * np.log10(lik / max_lik + 1e-300)) for gt, lik in likelihoods.items()}
    
    return likelihoods, posteriors, pl

print("=== Example: AD=20,22 DP=42 (clean heterozygous call) ===")
lik, post, pl = genotype_likelihood(ref_reads=20, alt_reads=22)
for gt in ['0/0','0/1','1/1']:
    print(f"  {gt}: likelihood={lik[gt]:.2e}  posterior={post[gt]:.4f}  PL={pl[gt]}")
print(f"  --> Called genotype: {max(post, key=post.get)}")

print("\n=== Example: AD=38,4 DP=42 (likely sequencing error, not real het) ===")
lik, post, pl = genotype_likelihood(ref_reads=38, alt_reads=4)
for gt in ['0/0','0/1','1/1']:
    print(f"  {gt}: likelihood={lik[gt]:.2e}  posterior={post[gt]:.4f}  PL={pl[gt]}")
print(f"  --> Called genotype: {max(post, key=post.get)}")
# Visualize genotype likelihood surfaces across different AD ratios
library(ggplot2)
alt_fracs <- seq(0, 1, 0.01)
dp <- 40

lik_00 <- dbinom(round(alt_fracs*dp), dp, 0.01)
lik_01 <- dbinom(round(alt_fracs*dp), dp, 0.50)
lik_11 <- dbinom(round(alt_fracs*dp), dp, 0.99)

gt_df <- data.frame(
  alt_frac = rep(alt_fracs, 3),
  likelihood = c(lik_00, lik_01, lik_11),
  Genotype = rep(c("0/0 (hom-ref)","0/1 (het)","1/1 (hom-alt)"), each=length(alt_fracs))
)

ggplot(gt_df, aes(x=alt_frac, y=likelihood, color=Genotype)) +
  geom_line(linewidth=1.3) +
  geom_vline(xintercept=c(0.05, 0.5, 0.95), linetype="dotted", color="grey50") +
  scale_color_manual(values=c("0/0 (hom-ref)"="#3498db","0/1 (het)"="#f39c12","1/1 (hom-alt)"="#e74c3c")) +
  labs(title="Genotype Likelihood as a Function of Observed Alt Allele Fraction",
       subtitle="DP=40; peaks at the fraction each genotype predicts",
       x="Observed Alt Allele Fraction", y="Likelihood P(D|G)") +
  theme_minimal(base_size=12)

Reference: Li H. (2011). “A statistical framework for SNP calling…” Bioinformatics 27(21):2987–2993 (defines the Bayesian genotype likelihood model used by GATK/samtools/bcftools). [Link]

18.5.2 16.3.2 Principal Component Analysis and Singular Value Decomposition — from Linear Algebra to UMAP

Pre-doctoral depth: PCA (used throughout Part 6 for scRNA-seq dimensionality reduction) is, at its core, an eigendecomposition of the covariance matrix, computed via Singular Value Decomposition (SVD): \[X = U \Sigma V^T\] where X is the (cells × genes) expression matrix, U’s columns are the principal component scores, and \(\Sigma\)’s diagonal contains the singular values (proportional to the variance explained by each component). Computing full SVD is O(min(n,p)² × max(n,p)); for scRNA-seq with p~20,000 genes and n~10,000 cells, truncated/randomized SVD (Halko et al. 2011) reduces this to near-linear time for the top k components — this is what scanpy.pp.pca and Seurat::RunPCA actually call under the hood (ARPACK or randomized solvers).

import numpy as np
from sklearn.decomposition import PCA, TruncatedSVD
import matplotlib.pyplot as plt

# Simulate a simplified single-cell expression matrix: 3 "cell types" with distinct expression programs
np.random.seed(42)
n_cells_per_type = 100
n_genes = 500

def simulate_celltype(mean_expr, n_cells, n_genes, noise=0.3):
    return np.random.lognormal(mean=mean_expr, sigma=noise, size=(n_cells, n_genes))

# 3 cell types with different mean expression programs (simplified — real data has gene-specific programs)
type_A = simulate_celltype(2.0, n_cells_per_type, n_genes)
type_B = simulate_celltype(2.0, n_cells_per_type, n_genes)
type_B[:, :50] *= 5   # cell-type-B-specific marker genes upregulated
type_C = simulate_celltype(2.0, n_cells_per_type, n_genes)
type_C[:, 50:100] *= 5  # cell-type-C-specific markers

X = np.vstack([type_A, type_B, type_C])
labels = ['A']*n_cells_per_type + ['B']*n_cells_per_type + ['C']*n_cells_per_type

# Log-normalize (standard scRNA-seq preprocessing)
X_norm = np.log1p(X)
X_scaled = (X_norm - X_norm.mean(axis=0)) / (X_norm.std(axis=0) + 1e-8)

# PCA via truncated SVD (what scanpy/Seurat use for speed on sparse matrices)
pca = PCA(n_components=10, svd_solver='randomized', random_state=42)
X_pca = pca.fit_transform(X_scaled)

print("Explained variance ratio (first 5 PCs):")
for i, var in enumerate(pca.explained_variance_ratio_[:5]):
    print(f"  PC{i+1}: {var*100:.2f}%")

print(f"\nCumulative variance explained (10 PCs): {pca.explained_variance_ratio_.sum()*100:.1f}%")

# Visualize PC1 vs PC2
plt.figure(figsize=(7,5))
colors = {'A':'#3498db','B':'#e74c3c','C':'#27ae60'}
for ct in ['A','B','C']:
    mask = [l==ct for l in labels]
    plt.scatter(X_pca[mask,0], X_pca[mask,1], label=f'Cell type {ct}', alpha=0.6, color=colors[ct])
plt.xlabel(f'PC1 ({pca.explained_variance_ratio_[0]*100:.1f}%)')
plt.ylabel(f'PC2 ({pca.explained_variance_ratio_[1]*100:.1f}%)')
plt.title('PCA of Simulated Single-Cell Expression (3 cell types)')
plt.legend()
plt.tight_layout()
plt.savefig('figures/pca_demo.png', dpi=150)
print("\nSaved PCA plot to figures/pca_demo.png")

Reference: Halko N., Martinsson P.G., Tropp J.A. (2011). “Finding structure with randomness: probabilistic algorithms for constructing approximate matrix decompositions.” SIAM Review 53(2):217–288 (randomized SVD used in scanpy/Seurat). [Link]

18.5.3 16.3.3 From Classical ML to Deep Learning: DeepVariant and AlphaFold Architecture

ml_arch <- data.frame(
  Model = c("Random Forest (VQSR-adjacent)","CNN (DeepVariant)","LSTM/RNN (early splicing models)",
             "Transformer (SpliceAI successor, Enformer)","Graph Neural Network (protein structure)",
             "Diffusion Model (RFdiffusion)","AlphaFold2 (Evoformer + Structure Module)"),
  Core_Mechanism = c(
    "Ensemble of decision trees; bagging + random feature subsets",
    "Convolutional filters slide over pileup 'image'; hierarchical feature learning",
    "Sequential hidden state updates; captures long-range dependence (superseded)",
    "Self-attention over sequence positions; captures very long-range regulatory context",
    "Message passing between graph nodes (residues); learns 3D spatial relationships",
    "Iterative denoising from Gaussian noise toward valid structures",
    "Multiple sequence alignment (MSA) + pairwise representation + iterative refinement"
  ),
  Genomics_Application = c("Variant filtering (older GATK)","Germline variant calling",
                             "Early splice-site prediction","Long-range gene regulation prediction (Enformer)",
                             "Protein structure prediction, drug binding","De novo protein design",
                             "Protein structure prediction from sequence"),
  Reference = c("Breiman 2001","Poplin 2018 NBT","Xiong 2015 Science",
                "Avsec 2021 Nat Methods","Jumper 2021 Nature (AF2 module)",
                "Watson 2023 Nature","Jumper 2021 Nature"),
  stringsAsFactors=FALSE
)
kable(ml_arch, caption="Machine learning architectures used across modern genomics") %>%
  kable_styling(bootstrap_options=c("striped","hover","condensed"), full_width=TRUE, font_size=11)
Machine learning architectures used across modern genomics
Model Core_Mechanism Genomics_Application Reference
Random Forest (VQSR-adjacent) Ensemble of decision trees; bagging + random feature subsets Variant filtering (older GATK) Breiman 2001
CNN (DeepVariant) Convolutional filters slide over pileup ‘image’; hierarchical feature learning Germline variant calling Poplin 2018 NBT
LSTM/RNN (early splicing models) Sequential hidden state updates; captures long-range dependence (superseded) Early splice-site prediction Xiong 2015 Science
Transformer (SpliceAI successor, Enformer) Self-attention over sequence positions; captures very long-range regulatory context Long-range gene regulation prediction (Enformer) Avsec 2021 Nat Methods
Graph Neural Network (protein structure) Message passing between graph nodes (residues); learns 3D spatial relationships Protein structure prediction, drug binding Jumper 2021 Nature (AF2 module)
Diffusion Model (RFdiffusion) Iterative denoising from Gaussian noise toward valid structures De novo protein design Watson 2023 Nature
AlphaFold2 (Evoformer + Structure Module) Multiple sequence alignment (MSA) + pairwise representation + iterative refinement Protein structure prediction from sequence Jumper 2021 Nature 
# Educational: minimal CNN architecture resembling DeepVariant's pileup classifier
# (Using pseudocode-level PyTorch — illustrates the architecture, not production-ready)

architecture_summary = """
DeepVariant CNN Architecture (Inception-v3 based, simplified view):

Input: RGB "pileup image" tensor, shape (channels=6, height=100 reads, width=221 positions)
  Channel 1: base identity (A/C/G/T encoded as intensity)
  Channel 2: base quality (Phred score, scaled 0-254)
  Channel 3: mapping quality
  Channel 4: strand (+ or -)
  Channel 5: supports variant vs. supports reference
  Channel 6: read is a duplicate

  |
  v
[Conv2D 3x3, 32 filters] -> [BatchNorm] -> [ReLU]
  |
  v
[Inception modules x N]  <- parallel conv branches (1x1, 3x3, 5x5) concatenated
  |                          this multi-scale view captures both local (single-base)
  v                          and regional (indel context) signal simultaneously
[Global Average Pooling]
  |
  v
[Fully Connected -> Softmax]
  |
  v
Output: P(0/0), P(0/1), P(1/1)  <- 3-class genotype probability
"""
print(architecture_summary)
import torch
import torch.nn as nn

class SimplifiedVariantCNN(nn.Module):
    """
    Simplified educational CNN inspired by DeepVariant's pileup classifier.
    Real DeepVariant uses a full Inception-v3; this demonstrates the core idea
    in <30 lines for pedagogical purposes.
    """
    def __init__(self, in_channels=6, n_classes=3):
        super().__init__()
        self.conv1 = nn.Conv2d(in_channels, 32, kernel_size=3, padding=1)
        self.bn1   = nn.BatchNorm2d(32)
        self.conv2 = nn.Conv2d(32, 64, kernel_size=3, padding=1)
        self.bn2   = nn.BatchNorm2d(64)
        self.pool  = nn.AdaptiveAvgPool2d((1,1))
        self.fc    = nn.Linear(64, n_classes)
        self.relu  = nn.ReLU()

    def forward(self, x):
        x = self.relu(self.bn1(self.conv1(x)))
        x = self.pool(self.relu(self.bn2(self.conv2(x))))
        x = x.view(x.size(0), -1)
        logits = self.fc(x)
        return torch.softmax(logits, dim=1)

# Demonstrate with a random "pileup image" batch
model = SimplifiedVariantCNN()
dummy_pileup = torch.randn(4, 6, 100, 221)   # batch=4, channels=6, reads=100, positions=221
genotype_probs = model(dummy_pileup)
print("Output shape:", genotype_probs.shape)  # (4, 3) -> P(0/0), P(0/1), P(1/1) per sample
print("Example probabilities (sample 0):", genotype_probs[0].detach().numpy().round(3))
print(f"\nModel parameters: {sum(p.numel() for p in model.parameters()):,}")
print("(Real DeepVariant Inception-v3 model has ~25M parameters)")

References: Poplin R. et al. (2018). Nature Biotechnology 36:983–987; Szegedy C. et al. (2016). “Rethinking the Inception Architecture.” CVPR (Inception-v3 backbone); Jumper J. et al. (2021). “Highly accurate protein structure prediction with AlphaFold.” Nature 596:583–589. [AlphaFold2]

18.5.4 16.3.4 Population Genetics — Hardy-Weinberg, F_ST, and PCA-based Ancestry

# Hardy-Weinberg Equilibrium testing
# p^2 + 2pq + q^2 = 1, where p = ref allele freq, q = alt allele freq
hwe_test <- function(n_AA, n_Aa, n_aa) {
  n_total <- n_AA + n_Aa + n_aa
  p <- (2*n_AA + n_Aa) / (2*n_total)   # ref allele frequency
  q <- 1 - p                            # alt allele frequency
  
  expected_AA <- p^2 * n_total
  expected_Aa <- 2*p*q * n_total
  expected_aa <- q^2 * n_total
  
  chi_sq <- sum((c(n_AA,n_Aa,n_aa) - c(expected_AA,expected_Aa,expected_aa))^2 / 
                c(expected_AA,expected_Aa,expected_aa))
  p_value <- 1 - pchisq(chi_sq, df=1)
  
  cat(sprintf("Observed:  AA=%d  Aa=%d  aa=%d\n", n_AA, n_Aa, n_aa))
  cat(sprintf("Expected:  AA=%.1f  Aa=%.1f  aa=%.1f\n", expected_AA, expected_Aa, expected_aa))
  cat(sprintf("Allele freq: p=%.3f  q=%.3f\n", p, q))
  cat(sprintf("Chi-sq = %.3f, p-value = %.4f\n", chi_sq, p_value))
  cat(sprintf("HWE %s (alpha=0.05)\n", ifelse(p_value > 0.05, "NOT rejected — consistent with HWE", 
                                                 "REJECTED — deviation from HWE (genotyping error? selection? population structure?)")))
  invisible(list(chi_sq=chi_sq, p_value=p_value))
}

cat("=== Example: variant consistent with HWE ===\n")
#> === Example: variant consistent with HWE ===
hwe_test(n_AA=810, n_Aa=180, n_aa=10)
#> Observed:  AA=810  Aa=180  aa=10
#> Expected:  AA=810.0  Aa=180.0  aa=10.0
#> Allele freq: p=0.900  q=0.100
#> Chi-sq = 0.000, p-value = 1.0000
#> HWE NOT rejected — consistent with HWE (alpha=0.05)
cat("\n=== Example: variant deviating from HWE (possible genotyping artifact) ===\n")
#> 
#> === Example: variant deviating from HWE (possible genotyping artifact) ===
hwe_test(n_AA=700, n_Aa=100, n_aa=200)
#> Observed:  AA=700  Aa=100  aa=200
#> Expected:  AA=562.5  Aa=375.0  aa=62.5
#> Allele freq: p=0.750  q=0.250
#> Chi-sq = 537.778, p-value = 0.0000
#> HWE REJECTED — deviation from HWE (genotyping error? selection? population structure?) (alpha=0.05)

Why HWE testing matters in QC: a variant that significantly deviates from Hardy-Weinberg equilibrium in a large, randomly-mating population sample is often flagged as a likely genotyping artifact (e.g., a paralogous region misaligned as a SNP) rather than true biology — a standard variant QC filter (--hwe 1e-6 in PLINK/vcftools) in GWAS pipelines.

# F_ST: population differentiation statistic
import numpy as np

def compute_fst(pop1_allele_freqs, pop2_allele_freqs, pop1_n, pop2_n):
    """
    Weir & Cockerham (1984) F_ST estimator (simplified single-locus version).
    F_ST measures genetic differentiation between populations:
      0 = no differentiation, 1 = complete differentiation (fixed different alleles)
    """
    p1, p2 = pop1_allele_freqs, pop2_allele_freqs
    n1, n2 = pop1_n, pop2_n
    n_total = n1 + n2
    
    p_bar = (n1*p1 + n2*p2) / n_total  # weighted mean allele frequency
    
    # Between-population variance
    s_squared = (n1*(p1-p_bar)**2 + n2*(p2-p_bar)**2) / n_total
    
    # Expected heterozygosity
    h_bar = 2 * p_bar * (1 - p_bar)
    
    fst = s_squared / (p_bar * (1 - p_bar)) if p_bar*(1-p_bar) > 0 else 0
    return fst

# Example: a variant with very different frequencies across populations (ancestry-informative marker)
fst_high = compute_fst(pop1_allele_freqs=0.05, pop2_allele_freqs=0.65, pop1_n=500, pop2_n=500)
print(f"F_ST (ancestry-informative marker): {fst_high:.4f}  <- high differentiation")

fst_low = compute_fst(pop1_allele_freqs=0.30, pop2_allele_freqs=0.32, pop1_n=500, pop2_n=500)
print(f"F_ST (typical neutral variant):     {fst_low:.4f}  <- low differentiation")

References: Weir B.S., Cockerham C.C. (1984). “Estimating F-statistics for the analysis of population structure.” Evolution 38(6):1358–1370. Price A.L. et al. (2006). “Principal components analysis corrects for stratification in genome-wide association studies.” Nature Genetics 38:904–909 (PCA for ancestry — the standard GWAS QC step). [Link]

18.6 16.4 TIER 4 — Systems, Software Engineering, and Reproducibility (Pre-Doctoral / Industry)

18.6.1 16.4.1 Complexity and Scale: What “Big Data” Means in Genomics

scale_table <- data.frame(
  Data_Type = c("Single WGS FASTQ pair (30x)","Single WGS BAM","Single WES BAM",
                 "1000 Genomes Project (full)","UK Biobank WGS (~500K genomes)",
                 "gnomAD v4 joint VCF","Single 10x scRNA-seq run (10K cells)",
                 "TCGA pan-cancer (multi-omic)","All of Us Research Program (target)"),
  Approx_Size = c("~90 GB","~90 GB (compressed CRAM: ~30GB)","~15 GB",
                   "~30 TB (raw)","~15-20 PB (raw)","~500 GB (VCF), ~1 PB source data",
                   "~5-10 GB (raw), ~500MB (matrix)","~2.5 PB (all data types)","~1 EB (projected)"),
  Compute_Note = c("~4-8 CPU-hours to align+call","—","~1-2 CPU-hours",
                    "Cloud/HPC required","Requires distributed cloud infrastructure (AWS/GCP)",
                    "Joint genotyping needs GenomicsDB sharding across intervals",
                    "GPU acceleration common for downstream ML",
                    "Requires federated/tiered storage strategy","Requires exabyte-scale federated architecture"),
  stringsAsFactors=FALSE
)
kable(scale_table, caption="Data scale in modern genomics — why distributed systems matter") %>%
  kable_styling(bootstrap_options=c("striped","hover","condensed"), full_width=TRUE, font_size=11)
Data scale in modern genomics — why distributed systems matter
Data_Type Approx_Size Compute_Note
Single WGS FASTQ pair (30x) ~90 GB ~4-8 CPU-hours to align+call
Single WGS BAM ~90 GB (compressed CRAM: ~30GB) —
Single WES BAM ~15 GB ~1-2 CPU-hours
1000 Genomes Project (full) ~30 TB (raw) Cloud/HPC required
UK Biobank WGS (~500K genomes) ~15-20 PB (raw) Requires distributed cloud infrastructure (AWS/GCP)
gnomAD v4 joint VCF ~500 GB (VCF), ~1 PB source data Joint genotyping needs GenomicsDB sharding across intervals
Single 10x scRNA-seq run (10K cells) ~5-10 GB (raw), ~500MB (matrix) GPU acceleration common for downstream ML
TCGA pan-cancer (multi-omic) ~2.5 PB (all data types) Requires federated/tiered storage strategy
All of Us Research Program (target) ~1 EB (projected) Requires exabyte-scale federated architecture

18.6.2 16.4.2 Distributed Computing Patterns in Genomics Pipelines

MapReduce-style scatter-gather is the dominant parallelization pattern in genomics (used by GATK’s -L intervals, Spark-based tools, and Nextflow’s channel model):

  1. Scatter: split the genome into N intervals (e.g., per-chromosome or fixed-size windows)
  2. Map: run the same operation (e.g., HaplotypeCaller) independently and in parallel on each interval — embarrassingly parallel, no inter-task communication needed
  3. Gather: concatenate/merge the N output files back into one genome-wide result (GatherVcfs, MergeSamFiles)

This pattern scales near-linearly with available compute nodes because genomic intervals are (mostly) independent — the main exception is structural variant calling near interval boundaries, handled via interval padding.

# Educational: scatter-gather pattern implementation (what Nextflow automates)
import concurrent.futures
import time

def process_genomic_interval(interval: tuple) -> dict:
    """Simulates variant calling on one genomic interval (e.g., HaplotypeCaller -L chr1:1-1000000)."""
    chrom, start, end = interval
    time.sleep(0.1)  # simulate compute time
    n_variants = (end - start) // 10000  # simulate variant density
    return {'interval': f"{chrom}:{start}-{end}", 'n_variants': n_variants}

def scatter_gather_pipeline(genome_length=250_000_000, n_chunks=25, max_workers=8):
    """
    Scatter: split genome into chunks
    Map: process each chunk in parallel
    Gather: aggregate results
    """
    chunk_size = genome_length // n_chunks
    intervals = [("chr1", i*chunk_size, (i+1)*chunk_size) for i in range(n_chunks)]
    
    print(f"Scattering into {n_chunks} intervals, processing with {max_workers} parallel workers...")
    start_time = time.time()
    
    # Map (parallel execution — this is what Nextflow's executor does across nodes)
    with concurrent.futures.ThreadPoolExecutor(max_workers=max_workers) as executor:
        results = list(executor.map(process_genomic_interval, intervals))
    
    elapsed = time.time() - start_time
    
    # Gather
    total_variants = sum(r['n_variants'] for r in results)
    
    print(f"Completed in {elapsed:.2f}s (vs. {n_chunks*0.1:.2f}s if run serially)")
    print(f"Speedup: {(n_chunks*0.1)/elapsed:.1f}x")
    print(f"Total variants (gathered): {total_variants:,}")
    return results

results = scatter_gather_pipeline()

18.6.3 16.4.3 Software Engineering for Reproducibility

Interview Q (any senior bioinformatics role): “How do you ensure a pipeline run in 2024 produces bit-identical results if re-run in 2028?”

Model Answer, in layers:
1. Containerization — pin every tool to an exact version inside Docker/Singularity images (a bwa-mem2:2.2.1 tag, not latest). Container digests (SHA256) can be pinned for absolute reproducibility.
2. Workflow provenance — Nextflow/Snakemake record every parameter, container, and input file hash used per run (pipeline_info/ in nf-core).
3. Reference data versioning — pin exact reference genome build + annotation GTF version (GENCODE version numbers matter — gene models change between releases).
4. Deterministic algorithms — some tools (e.g., multi-threaded assemblers) are not fully deterministic across thread counts; document and pin thread counts if bit-identical output is required, or use tools with guaranteed determinism.
5. Version control — pipeline code itself in git, tagged releases, not “the version on my laptop.”
6. Environment capture — conda env export, renv::snapshot() (R), pip freeze / poetry.lock (Python) for any custom scripts outside containers.

# Example: fully pinned, reproducible pipeline invocation
nextflow run nf-core/sarek \
  -revision 3.4.0 \                                    # exact pipeline version (git tag)
  -profile docker \
  -with-docker \
  --input samplesheet.csv \
  --genome GATK.GRCh38 \                               # pinned reference build
  --igenomes_base s3://ngi-igenomes/igenomes/ \        # pinned reference source
  --tools haplotypecaller,vep \
  --vep_cache_version 110 \                             # pinned annotation version
  --outdir results/ \
  -with-report execution_report_$(date +%Y%m%d).html \
  -with-trace execution_trace_$(date +%Y%m%d).txt

# Container digest pinning (absolute reproducibility, not just tag)
# docker pull broadinstitute/gatk@sha256:a1b2c3...
cat("
# R reproducibility: renv (equivalent to Python's venv/poetry)
renv::init()                  # create isolated project library
renv::snapshot()              # lock exact package versions -> renv.lock
renv::restore()                # reinstall exact versions on a new machine

# Session info should always be included in any analysis report:
")
#> 
#> # R reproducibility: renv (equivalent to Python's venv/poetry)
#> renv::init()                  # create isolated project library
#> renv::snapshot()              # lock exact package versions -> renv.lock
#> renv::restore()                # reinstall exact versions on a new machine
#> 
#> # Session info should always be included in any analysis report:
sessionInfo()
#> R version 4.6.1 (2026-06-24 ucrt)
#> Platform: x86_64-w64-mingw32/x64
#> Running under: Windows 11 x64 (build 26200)
#> 
#> Matrix products: default
#>   LAPACK version 3.12.1
#> 
#> locale:
#> [1] LC_COLLATE=English_India.utf8  LC_CTYPE=English_India.utf8   
#> [3] LC_MONETARY=English_India.utf8 LC_NUMERIC=C                  
#> [5] LC_TIME=English_India.utf8    
#> 
#> time zone: Asia/Calcutta
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats4    stats     graphics  grDevices utils     datasets  methods  
#> [8] base     
#> 
#> other attached packages:
#>  [1] Biostrings_2.80.1           XVector_0.52.0             
#>  [3] Seurat_5.5.1                SeuratObject_5.4.0         
#>  [5] sp_2.2-3                    org.Hs.eg.db_3.23.1        
#>  [7] AnnotationDbi_1.74.0        clusterProfiler_4.20.0     
#>  [9] DESeq2_1.52.0               SummarizedExperiment_1.42.0
#> [11] Biobase_2.72.0              MatrixGenerics_1.24.0      
#> [13] matrixStats_1.5.0           GenomicRanges_1.64.0       
#> [15] Seqinfo_1.2.0               IRanges_2.46.0             
#> [17] S4Vectors_0.50.1            BiocGenerics_0.58.1        
#> [19] generics_0.1.4              apeglm_1.34.0              
#> [21] ggrepel_0.9.8               gridExtra_2.3.1            
#> [23] RColorBrewer_1.1-3          kableExtra_1.4.1           
#> [25] knitr_1.51                  scales_1.4.0               
#> [27] tidyr_1.3.2                 dplyr_1.2.1                
#> [29] ggplot2_4.0.3              
#> 
#> loaded via a namespace (and not attached):
#>   [1] fs_2.1.0                spatstat.sparse_3.2-0   enrichplot_1.32.0      
#>   [4] httr_1.4.8              numDeriv_2016.8-1.1     tools_4.6.1            
#>   [7] sctransform_0.4.3       R6_2.6.1                uwot_0.2.4             
#>  [10] lazyeval_0.2.3          mgcv_1.9-4              withr_3.0.3            
#>  [13] progressr_1.0.0         cli_3.6.6               textshaping_1.0.5      
#>  [16] spatstat.explore_3.8-2  fastDummies_1.7.6       scatterpie_0.2.6       
#>  [19] labeling_0.4.3          sass_0.4.10             mvtnorm_1.4-2          
#>  [22] S7_0.2.2                spatstat.data_3.1-9     ggridges_0.5.7         
#>  [25] pbapply_1.7-4           systemfonts_1.3.2       yulab.utils_0.2.4      
#>  [28] gson_0.2.1              DOSE_4.6.0              svglite_2.2.2          
#>  [31] parallelly_1.48.0       bbmle_1.0.25.1          rstudioapi_0.19.0      
#>  [34] RSQLite_3.53.3          gridGraphics_0.5-1      ica_1.0-3              
#>  [37] spatstat.random_3.5-1   GO.db_3.23.1            Matrix_1.7-5           
#>  [40] abind_1.4-8             lifecycle_1.0.5         yaml_2.3.12            
#>  [43] qvalue_2.44.0           SparseArray_1.12.2      Rtsne_0.17             
#>  [46] grid_4.6.1              blob_1.3.0              promises_1.5.0         
#>  [49] crayon_1.5.3            bdsmatrix_1.3-7         miniUI_0.1.2           
#>  [52] ggtangle_0.1.2          lattice_0.22-9          cowplot_1.2.0          
#>  [55] KEGGREST_1.52.2         pillar_1.11.1           future.apply_1.20.2    
#>  [58] codetools_0.2-20        glue_1.8.1              ggiraph_0.9.6          
#>  [61] ggfun_0.2.1             spatstat.univar_3.2-0   fontLiberation_0.1.0   
#>  [64] data.table_1.18.4       vctrs_0.7.3             png_0.1-9              
#>  [67] treeio_1.36.1           spam_2.11-4             gtable_0.3.6           
#>  [70] emdbook_1.3.14          cachem_1.1.0            xfun_0.60              
#>  [73] S4Arrays_1.12.0         mime_0.13               coda_0.19-4.1          
#>  [76] survival_3.8-6          aisdk_1.4.12            fitdistrplus_1.2-6     
#>  [79] ROCR_1.0-12             nlme_3.1-169            ggtree_4.2.0           
#>  [82] bit64_4.8.4             fontquiver_0.2.1        RcppAnnoy_0.0.23       
#>  [85] bslib_0.12.0            irlba_2.3.7             KernSmooth_2.23-26     
#>  [88] otel_0.2.0              DBI_1.3.0               tidyselect_1.2.1       
#>  [91] processx_3.9.0          bit_4.6.0               compiler_4.6.1         
#>  [94] httr2_1.3.0             xml2_1.6.0              fontBitstreamVera_0.1.1
#>  [97] DelayedArray_0.38.2     plotly_4.12.1           lmtest_0.9-40          
#> [100] callr_3.8.0             rappdirs_0.3.4          goftest_1.2-3          
#> [103] stringr_1.6.0           digest_0.6.39           spatstat.utils_3.2-4   
#> [106] rmarkdown_2.31          htmltools_0.5.9         pkgconfig_2.0.3        
#> [109] fastmap_1.2.0           rlang_1.3.0             htmlwidgets_1.6.4      
#> [112] shiny_1.14.0            farver_2.1.2            jquerylib_0.1.4        
#> [115] zoo_1.9-0               jsonlite_2.0.0          BiocParallel_1.46.0    
#> [118] GOSemSim_2.38.3         magrittr_2.0.5          ggplotify_0.1.3        
#> [121] dotCall64_1.2           patchwork_1.3.2         Rcpp_1.1.2             
#> [124] ape_5.8-1               ggnewscale_0.5.2        gdtools_0.5.1          
#> [127] reticulate_1.46.0       stringi_1.8.9           MASS_7.3-65            
#> [130] plyr_1.8.9              parallel_4.6.1          listenv_1.0.0          
#> [133] deldir_2.0-4            splines_4.6.1           tensor_1.5.1           
#> [136] locfit_1.5-9.12         igraph_2.3.3            spatstat.geom_3.8-2    
#> [139] enrichit_0.2.1          RcppHNSW_0.7.0          reshape2_1.4.5         
#> [142] evaluate_1.0.5          tweenr_2.0.3            httpuv_1.6.17          
#> [145] RANN_2.6.2              purrr_1.2.2             polyclip_1.10-7        
#> [148] future_1.75.0           scattermore_1.2         ggforce_0.5.0          
#> [151] xtable_1.8-8            RSpectra_0.16-2         tidytree_0.4.8         
#> [154] tidydr_0.0.6            later_1.4.8             viridisLite_0.4.3      
#> [157] tibble_3.3.1            aplot_0.3.1             memoise_2.0.1          
#> [160] cluster_2.1.8.2         globals_0.19.1

18.6.4 16.4.4 Testing Bioinformatics Code — Unit Tests for Genomics Functions

import unittest

def compute_gc_content(seq: str) -> float:
    """Compute GC content as a fraction."""
    if not seq: return 0.0
    seq = seq.upper()
    gc = seq.count('G') + seq.count('C')
    return gc / len(seq)

def parse_vcf_genotype(gt_field: str) -> tuple:
    """Parse a VCF GT field like '0/1' or '1|1' into a tuple of alleles."""
    sep = '|' if '|' in gt_field else '/'
    return tuple(int(a) if a != '.' else None for a in gt_field.split(sep))

class TestGenomicsFunctions(unittest.TestCase):
    """
    Unit tests for core genomics utility functions.
    In production: pytest + hypothesis (property-based testing) for edge cases.
    """
    
    def test_gc_content_basic(self):
        self.assertAlmostEqual(compute_gc_content("GCGC"), 1.0)
        self.assertAlmostEqual(compute_gc_content("ATAT"), 0.0)
        self.assertAlmostEqual(compute_gc_content("ATGC"), 0.5)
    
    def test_gc_content_empty(self):
        self.assertEqual(compute_gc_content(""), 0.0)
    
    def test_gc_content_lowercase(self):
        # Edge case: soft-masked (lowercase) repeat regions in reference genomes
        self.assertAlmostEqual(compute_gc_content("gcgc"), 1.0)
    
    def test_gt_parsing_unphased_het(self):
        self.assertEqual(parse_vcf_genotype("0/1"), (0, 1))
    
    def test_gt_parsing_phased(self):
        self.assertEqual(parse_vcf_genotype("1|0"), (1, 0))
    
    def test_gt_parsing_missing(self):
        # Edge case: missing genotype calls (low coverage)
        self.assertEqual(parse_vcf_genotype("./."), (None, None))
    
    def test_gt_parsing_hom_alt(self):
        self.assertEqual(parse_vcf_genotype("1/1"), (1, 1))

# Run tests
suite = unittest.TestLoader().loadTestsFromTestCase(TestGenomicsFunctions)
runner = unittest.TextTestRunner(verbosity=2)
runner.run(suite)

Industry note: production genomics pipelines (Illumina DRAGEN, GATK, nf-core modules) all maintain extensive test suites — nf-core requires nf-test coverage for every module (Part 12, Step 9) before it can be merged into the community pipeline library. Interviewers at Illumina, 10x Genomics, and DNAnexus commonly ask candidates to write unit tests alongside any coding exercise, not just working code.

18.7 16.5 TIER 5 — Pre-Doctoral Synthesis: Open Research Questions

A pre-doctoral bioinformatician should be able to discuss open problems, not just solved pipelines. Selected active research frontiers as of 2024–2026:

  1. Pangenome references — GRCh38/T2T-CHM13 are single linear references; the Human Pangenome Reference Consortium’s graph-based pangenome (47+ diverse haplotypes) aims to reduce reference bias, especially in structurally variable and historically underrepresented populations.
    Reference: Liao W.W. et al. (2023). “A draft human pangenome reference.” Nature 617:312–324. [Link]

  2. Single-cell foundation models — large pretrained transformer models (scGPT, Geneformer) trained on tens of millions of cells, aiming for zero-shot cell-type annotation, perturbation prediction, and gene-network inference, analogous to LLMs in NLP.
    Reference: Cui H. et al. (2024). “scGPT: toward building a foundation model for single-cell multi-omics using generative AI.” Nature Methods 21:1470–1480. [Link]

  3. Long-read clinical adoption — as Nanopore/PacBio accuracy approaches short-read levels, the field is shifting toward long-read-first clinical pipelines (resolving repeat expansions, complex SVs, and phasing in one assay) — an active area of pipeline redesign at major clinical labs as of 2024-2026.

  4. Variant effect prediction via protein language models — AlphaMissense and ESM (Evolutionary Scale Modeling) apply protein-language-model architectures (trained via masked-token prediction, analogous to BERT) to predict pathogenicity without labeled training data, purely from evolutionary sequence context.
    Reference: Cheng J. et al. (2023). “Accurate proteome-wide missense variant effect prediction with AlphaMissense.” Science 381:eadg7492. [Link]

  5. Spatial multi-omics integration — combining spatial transcriptomics, spatial proteomics, and histology into unified computational frameworks remains an unsolved integration challenge (batch effects across modalities are worse than within a single modality).

# Visualize the growth of genomics data and open research area investment (illustrative)
research_areas <- data.frame(
  Area = c("Pangenome\nReferences","Single-cell\nFoundation Models","Long-read\nClinical Pipelines",
           "Protein Language\nModels","Spatial Multi-omics\nIntegration","AI-driven\nDrug Discovery"),
  Maturity_2024 = c(3, 2, 4, 5, 2, 3),
  Projected_2028 = c(7, 7, 8, 8, 6, 7)
)

research_long <- research_areas %>%
  tidyr::pivot_longer(cols=c(Maturity_2024, Projected_2028), names_to="Year", values_to="Maturity")

ggplot(research_long, aes(x=Area, y=Maturity, fill=Year)) +
  geom_bar(stat="identity", position="dodge", width=0.6) +
  scale_fill_manual(values=c("Maturity_2024"="#95a5a6","Projected_2028"="#3498db"),
                    labels=c("2024 (current)","2028 (projected)")) +
  labs(title="Open Research Frontiers in Computational Genomics",
       subtitle="Illustrative maturity scale (1=nascent, 10=production-standard) — for orientation, not forecasting",
       x="", y="Field Maturity (illustrative)", fill="") +
  theme_minimal(base_size=11) +
  theme(axis.text.x=element_text(angle=20, hjust=1), legend.position="bottom")

18.8 16.6 Self-Assessment Question Bank (Beginner → Pre-Doctoral)

Beginner:
1. Why must a DNA aligner check both the forward and reverse-complement strand of every read?
2. What does Q30 mean, and why is it the de facto industry standard threshold?
3. Explain, in one sentence, why RNA-seq alignment requires a “splice-aware” aligner while DNA-seq does not.

Intermediate:
4. Derive why BWA’s FM-index search is O(m) independent of genome size, while naive search is O(nm).
5. Why does DESeq2 use a Negative Binomial model instead of Poisson, and what does “dispersion shrinkage” solve?
6. Explain how the Benjamini-Hochberg procedure differs mathematically from Bonferroni correction, and why genomics prefers BH.

Pre-Doctoral / Advanced:
7. Sketch how you would design a new statistical test for detecting differentially abundant cell types between two conditions in scRNA-seq data, accounting for the compositional (sum-to-one) nature of cell-type proportions.
8. Compare the computational complexity and biological assumptions of De Bruijn-graph assembly vs. overlap-layout-consensus (OLC) assembly. Under what read-length/error-rate regime does each dominate?
9. Given a pretrained protein language model (e.g., ESM2) and a VCF of missense variants, outline an end-to-end computational pipeline to prioritize variants for functional follow-up, including how you would benchmark it against ClinVar.
10. A collaborator proposes training a single foundation model across bulk RNA-seq, scRNA-seq, and spatial transcriptomics jointly. What are the key technical obstacles (batch effects, resolution mismatch, platform-specific noise models), and how might you address at least one?

19 PART 16 — Beginner → Pre-Doctoral Deep Dive for CS & Bioinformatics Graduates

This part is a self-contained, cumulative curriculum. It starts from “what is a base pair” and ends at graduate-level algorithmic and statistical machinery (suffix arrays, HMMs, EM, variational inference, GNNs for genomics) — written for someone with a computer science background who wants full mastery of the computational side of genomics, not just tool usage. Every subsection increases in depth; skip ahead if a level is already familiar, but each layer assumes the previous one.

19.1 16.0 Roadmap

levels_df <- data.frame(
  Level = factor(c("L0 Beginner","L1 Undergrad","L2 Advanced Undergrad",
                    "L3 MSc","L4 Pre-Doctoral"),
                 levels=c("L0 Beginner","L1 Undergrad","L2 Advanced Undergrad",
                          "L3 MSc","L4 Pre-Doctoral")),
  Topics = c("DNA/RNA/protein, central dogma, what a FASTQ is",
             "Strings, arrays, hash tables, O-notation, basic probability",
             "Suffix arrays/trees, DP alignment, HMMs, hypothesis testing",
             "GATK internals, EM algorithms, negative binomial GLMs, graph clustering",
             "Pair-HMMs, variational inference, GNNs, population genetics theory, causal inference"),
  X = 1:5, Y = rep(1,5)
)
ggplot(levels_df, aes(x=X, y=Y)) +
  geom_segment(aes(xend=X+0.001, yend=Y), linewidth=0) +
  geom_point(size=14, color="#2980b9") +
  geom_text(aes(label=1:5), color="white", fontface="bold", size=5) +
  geom_segment(data=data.frame(x=1:4, xend=2:5), aes(x=x,xend=xend,y=1,yend=1),
               inherit.aes=FALSE, arrow=arrow(length=unit(0.3,"cm")), color="grey50") +
  geom_text(aes(label=Level, y=1.15), size=3.5, fontface="bold") +
  geom_text(data=levels_df, aes(label=stringr::str_wrap(Topics,22), y=0.8),
            size=2.6, vjust=1, lineheight=0.85) +
  ylim(0.3,1.3) +
  theme_void() +
  labs(title="Curriculum Progression: Beginner → Pre-Doctoral Bioinformatics") +
  theme(plot.title=element_text(hjust=0.5, face="bold", size=13))
Learning progression from beginner to pre-doctoral level

Learning progression from beginner to pre-doctoral level

19.2 16.1 Level 0 — Absolute Beginner: Molecular Biology for Programmers

Historical anchor: Watson & Crick’s 1953 double-helix structure (built on Rosalind Franklin’s X-ray diffraction “Photo 51”) established that DNA is a self-complementary code: A pairs with T, G pairs with C, held by hydrogen bonds, arranged antiparallel (5′→3′ on one strand, 3′→5′ on the other).
Reference: Watson J.D., Crick F.H.C. (1953). “Molecular structure of nucleic acids.” Nature 171:737–738. [Link]

The Central Dogma (Crick, 1958): DNA → RNA → Protein.

dogma <- data.frame(x=c(1,2,3), y=c(1,1,1), label=c("DNA\n(4-letter code)","RNA\n(transcription)","Protein\n(translation)"))
ggplot(dogma, aes(x=x,y=y)) +
  geom_point(size=30, color=c("#3498db","#e67e22","#27ae60")) +
  geom_text(aes(label=label), color="white", fontface="bold", size=3.2) +
  geom_segment(data=data.frame(x=c(1,2),xend=c(2,3),y=c(1,1),yend=c(1,1)),
               aes(x=x,xend=xend,y=y,yend=yend), inherit.aes=FALSE,
               arrow=arrow(length=unit(0.4,"cm")), linewidth=1.2) +
  annotate("text", x=1.5, y=1.35, label="Transcription\n(RNA Pol II)", size=3) +
  annotate("text", x=2.5, y=1.35, label="Translation\n(Ribosome)", size=3) +
  xlim(0.5,3.5) + ylim(0.5,1.6) +
  theme_void() +
  labs(title="The Central Dogma of Molecular Biology (Crick, 1958)")

Think of it as a computer scientist:

Biology CS analogy
DNA (genome) Source code stored on disk (immutable master copy)
Gene A function definition within the source
RNA (transcript) Compiled bytecode / intermediate representation
Ribosome/translation The interpreter/runtime executing bytecode
Protein The running program’s actual behavior/output
Mutation A single-character edit to source code — may be a no-op, a warning, or a crash
Epigenetics (methylation) Environment variables/config flags that change which functions get called, without changing the code itself

The alphabet: DNA uses 4 symbols (A,C,G,T); RNA substitutes U for T. Each amino acid is coded by a codon — a 3-letter “word” from the 4-letter alphabet, giving 4³=64 possible codons for only 20 amino acids (degenerate code — most amino acids have multiple synonymous codons, which is why synonymous mutations are usually silent).

# The genetic code as a Python dictionary — first thing every bioinformatician implements once
CODON_TABLE = {
    'TTT':'F','TTC':'F','TTA':'L','TTG':'L','CTT':'L','CTC':'L','CTA':'L','CTG':'L',
    'ATT':'I','ATC':'I','ATA':'I','ATG':'M','GTT':'V','GTC':'V','GTA':'V','GTG':'V',
    'TCT':'S','TCC':'S','TCA':'S','TCG':'S','CCT':'P','CCC':'P','CCA':'P','CCG':'P',
    'ACT':'T','ACC':'T','ACA':'T','ACG':'T','GCT':'A','GCC':'A','GCA':'A','GCG':'A',
    'TAT':'Y','TAC':'Y','TAA':'*','TAG':'*','CAT':'H','CAC':'H','CAA':'Q','CAG':'Q',
    'AAT':'N','AAC':'N','AAA':'K','AAG':'K','GAT':'D','GAC':'D','GAA':'E','GAG':'E',
    'TGT':'C','TGC':'C','TGA':'*','TGG':'W','CGT':'R','CGC':'R','CGA':'R','CGG':'R',
    'AGT':'S','AGC':'S','AGA':'R','AGG':'R','GGT':'G','GGC':'G','GGA':'G','GGG':'G'
}  # '*' = stop codon

def translate(rna_or_dna_seq):
    """Translate a nucleotide sequence into a protein sequence. O(n) time."""
    seq = rna_or_dna_seq.upper().replace('U', 'T')
    protein = []
    for i in range(0, len(seq) - 2, 3):
        codon = seq[i:i+3]
        aa = CODON_TABLE.get(codon, 'X')  # X = unknown/ambiguous
        if aa == '*':
            break  # stop codon terminates translation
        protein.append(aa)
    return ''.join(protein)

print(translate("ATGGCCATTGTAATGGGCCGCTGAAAG"))
# ATG(M) GCC(A) ATT(I) GTA(V) ATG(M) GGC(G) CGC(R) TGA(*stop)
# -> "MAIVMGR"

Why bioinformaticians must know reverse complementation:

def reverse_complement(seq):
    """
    DNA is double-stranded and antiparallel — a gene can be encoded
    on either strand. Aligners must consider both.
    Time: O(n), Space: O(n)
    """
    complement = {'A':'T','T':'A','C':'G','G':'C','N':'N'}
    return ''.join(complement[base] for base in reversed(seq.upper()))

seq = "ATGGCCATT"
print(f"Original:            5'-{seq}-3'")
print(f"Reverse complement:  5'-{reverse_complement(seq)}-3'")
# This is why aligners test both the forward AND reverse-complement strand for every read

19.3 16.2 Level 1 — Undergraduate CS Foundations Applied to Genomics

19.3.1 16.2.1 Strings, Arrays, and Why Genomics Is a String-Processing Field at Scale

A human genome is a string of 3.1 × 10⁹ characters over a 4-letter alphabet (ignoring N’s and lowercase soft-masking). Every “alignment,” “assembly,” or “search” problem in genomics reduces, at its core, to a classic CS string-processing problem — but at a scale where naive algorithms are computationally intractable.

scale_df <- data.frame(
  Object = c("Human genome","Human exome (protein-coding)","Single WGS FASTQ read",
              "One lane of NovaSeq X (WGS)","1000 Genomes cohort VCF","gnomAD v4 database"),
  Approx_Size = c("3.1 Gb (~3.1×10⁹ chars)","~30 Mb (~1% of genome)","150 bp",
                   "~1.5 Tb raw data","~200 GB (compressed)","~15 TB (all variant data)"),
  CS_Analogy = c("A single 3GB text file with no line breaks",
                  "A ~30MB file: the 'real code' hidden in the 3GB repo",
                  "One 150-character string — like a tweet",
                  "Roughly 500,000 human genomes' worth of raw text per day, per instrument",
                  "A sparse matrix with 2,504 columns (samples) x 84M rows (variants)",
                  "Bigger than the entire English Wikipedia text dump, many times over"),
  stringsAsFactors=FALSE
)
kable(scale_df, caption="Scale intuition: genomics data sizes vs. familiar CS objects") %>%
  kable_styling(bootstrap_options=c("striped","hover"), full_width=TRUE)
Scale intuition: genomics data sizes vs. familiar CS objects
Object Approx_Size CS_Analogy
Human genome 3.1 Gb (~3.1×10⁹ chars) A single 3GB text file with no line breaks
Human exome (protein-coding) ~30 Mb (~1% of genome) A ~30MB file: the ‘real code’ hidden in the 3GB repo
Single WGS FASTQ read 150 bp One 150-character string — like a tweet
One lane of NovaSeq X (WGS) ~1.5 Tb raw data Roughly 500,000 human genomes’ worth of raw text per day, per instrument
1000 Genomes cohort VCF ~200 GB (compressed) A sparse matrix with 2,504 columns (samples) x 84M rows (variants)
gnomAD v4 database ~15 TB (all variant data) Bigger than the entire English Wikipedia text dump, many times over

Complexity intuition every bioinformatics MSc/PhD student must internalize:

Algorithm Naive complexity Genomic-scale feasibility (n = 3×10⁹)
Naive substring search O(n·m) per query Infeasible (~10¹¹ ops/read × millions of reads)
Smith-Waterman (full DP) O(n·m) per pair Infeasible genome-wide; used only locally after seeding
Suffix array construction O(n log n) (or O(n) with SA-IS) Feasible once, reused for all queries
FM-index backward search O(m) per query, independent of n This is why BWA/Bowtie work
Hash table k-mer lookup O(m) average per query Basis of Salmon/kallisto/minimap2 seeding

19.3.2 16.2.2 Hash Tables in Genomics: k-mers

A k-mer is a substring of length k. Sequencing tools use hash tables mapping k-mer → genomic position(s) as an O(1)-average-lookup index.

def build_kmer_index(reference_seq, k=20):
    """
    Build a hash table: k-mer -> list of starting positions.
    Space: O(n) k-mers, each O(k) — total O(nk), but with hashing, effectively O(n).
    Time to build: O(n).
    This is the core data structure behind minimizer-based aligners (minimap2) 
    and k-mer counters (Jellyfish, KMC).
    """
    index = {}
    for i in range(len(reference_seq) - k + 1):
        kmer = reference_seq[i:i+k]
        index.setdefault(kmer, []).append(i)
    return index

def query_kmer(index, kmer):
    """O(1) average-case lookup — the entire point of using a hash table here."""
    return index.get(kmer, [])

reference = "ACGTACGTTGCATGCATGCA"
idx = build_kmer_index(reference, k=4)
print(f"Positions of 'ACGT': {query_kmer(idx, 'ACGT')}")
print(f"Positions of 'TGCA': {query_kmer(idx, 'TGCA')}")
print(f"Index size: {len(idx)} unique 4-mers out of {len(reference)-3} total windows")

Why k-mer size matters (a bias-variance tradeoff, in ML terms): - Small k (e.g., k=8): many false-positive hits (k-mer appears everywhere) → low specificity, but tolerant of sequencing errors - Large k (e.g., k=31, common in de Bruijn graph assemblers): highly specific, unique in most of genome, but a single sequencing error destroys the k-mer entirely (high variance to noise) - This exact tradeoff reappears in assembler design (SPAdes, Velvet — choice of k), in Kraken2 taxonomic classification, and in scRNA-seq pseudoalignment (Salmon/kallisto)

Reference: Marçais G., Kingsford C. (2011). “A fast, lock-free approach for efficient parallel counting of occurrences of k-mers.” Bioinformatics 27(6):764–770 (Jellyfish). [Link]

19.3.3 16.2.3 Dynamic Programming: Sequence Alignment From First Principles

This is the single most important CS concept in classical bioinformatics.

Needleman-Wunsch (global alignment, 1970) — the original application of dynamic programming to biology, predating its common use in CS algorithms courses.

import numpy as np

def needleman_wunsch(seq1, seq2, match=1, mismatch=-1, gap=-2):
    """
    Global alignment via dynamic programming.
    Time: O(n*m), Space: O(n*m) [can be reduced to O(min(n,m)) with Hirschberg's algorithm]
    
    Recurrence:
      F(i,j) = max(
          F(i-1,j-1) + score(seq1[i],seq2[j]),   # match/mismatch (diagonal)
          F(i-1,j)   + gap,                       # deletion (up)
          F(i,j-1)   + gap                        # insertion (left)
      )
    """
    n, m = len(seq1), len(seq2)
    F = np.zeros((n+1, m+1), dtype=int)
    
    # Initialize borders (cost of aligning to all-gaps)
    for i in range(n+1): F[i][0] = i * gap
    for j in range(m+1): F[0][j] = j * gap
    
    # Fill DP table
    for i in range(1, n+1):
        for j in range(1, m+1):
            score = match if seq1[i-1] == seq2[j-1] else mismatch
            F[i][j] = max(
                F[i-1][j-1] + score,   # diagonal
                F[i-1][j] + gap,       # up
                F[i][j-1] + gap        # left
            )
    
    # Traceback to recover the actual alignment
    align1, align2 = "", ""
    i, j = n, m
    while i > 0 or j > 0:
        if i > 0 and j > 0 and F[i][j] == F[i-1][j-1] + (match if seq1[i-1]==seq2[j-1] else mismatch):
            align1 = seq1[i-1] + align1
            align2 = seq2[j-1] + align2
            i -= 1; j -= 1
        elif i > 0 and F[i][j] == F[i-1][j] + gap:
            align1 = seq1[i-1] + align1
            align2 = "-" + align2
            i -= 1
        else:
            align1 = "-" + align1
            align2 = seq2[j-1] + align2
            j -= 1
    
    return F[n][m], align1, align2

score, a1, a2 = needleman_wunsch("GATTACA", "GCATGCU")
print(f"Alignment score: {score}")
print(f"Seq1: {a1}")
print(f"Seq2: {a2}")

Smith-Waterman (local alignment, 1981) — the key modification: clamp negative scores to 0, allowing the traceback to start anywhere, not just the bottom-right corner. This is the algorithm BWA-MEM uses to extend seeds found via the FM-index.

def smith_waterman(seq1, seq2, match=2, mismatch=-1, gap=-2):
    """
    Local alignment — finds the best-scoring SUBSEQUENCE alignment, not the whole strings.
    Key change from Needleman-Wunsch: F(i,j) = max(0, ...) — negative scores reset to zero,
    which allows alignments to 'restart' anywhere, ignoring non-matching flanks.
    This is exactly what you want when aligning a 150bp read to a 3Gb genome:
    you don't care about the whole genome matching, only the 150bp window that does.
    """
    n, m = len(seq1), len(seq2)
    F = [[0]*(m+1) for _ in range(n+1)]
    max_score, max_pos = 0, (0,0)
    
    for i in range(1, n+1):
        for j in range(1, m+1):
            score = match if seq1[i-1] == seq2[j-1] else mismatch
            F[i][j] = max(0,  # <-- the critical difference from global alignment
                          F[i-1][j-1] + score,
                          F[i-1][j] + gap,
                          F[i][j-1] + gap)
            if F[i][j] > max_score:
                max_score, max_pos = F[i][j], (i,j)
    
    return max_score, max_pos

score, pos = smith_waterman("AGCACACAGATCCCC", "ACACACTA")
print(f"Best local alignment score: {score} ending at reference position {pos}")

References: Needleman S.B., Wunsch C.D. (1970). “A general method applicable to the search for similarities in the amino acid sequence of two proteins.” JMB 48(3):443–453. Smith T.F., Waterman M.S. (1981). “Identification of common molecular subsequences.” JMB 147(1):195–197.

19.3.4 16.2.4 Basic Probability for NGS: Binomial and the Birthday Problem of Barcodes

# Lander-Waterman model: probability a genomic position has 0 coverage
# given N reads of length L covering a genome of size G
# Expected coverage lambda = N*L/G
# P(position uncovered) = e^(-lambda)  [Poisson approximation to Binomial]

G <- 3.1e9        # genome size
L <- 150          # read length
lambdas <- seq(1, 50, 1)

p_uncovered <- exp(-lambdas)
frac_genome_covered <- 1 - p_uncovered

lw_df <- data.frame(coverage=lambdas, p_gap=p_uncovered, frac_covered=frac_genome_covered)

ggplot(lw_df, aes(x=coverage, y=frac_covered*100)) +
  geom_line(color="#2980b9", linewidth=1.3) +
  geom_hline(yintercept=99, linetype="dashed", color="red") +
  geom_vline(xintercept=lw_df$coverage[which.min(abs(lw_df$frac_covered - 0.99))],
             linetype="dashed", color="red") +
  labs(title="Lander-Waterman Model: Genome Coverage vs. Sequencing Depth",
       subtitle="P(position covered) = 1 - e^(-lambda); explains why 30x is the WGS standard",
       x="Average Coverage (lambda = N*L/G)", y="% Genome Covered (Poisson model)") +
  theme_minimal(base_size=12)

cat(sprintf("At 30x coverage: %.4f%% of genome expected covered (Poisson model, ignoring bias)\n",
            (1-exp(-30))*100))
#> At 30x coverage: 100.0000% of genome expected covered (Poisson model, ignoring bias)
cat(sprintf("At 10x coverage: %.4f%% of genome expected covered\n", (1-exp(-10))*100))
#> At 10x coverage: 99.9955% of genome expected covered

Reference: Lander E.S., Waterman M.S. (1988). “Genomic mapping by fingerprinting random clones: a mathematical analysis.” Genomics 2(3):231–239 — the original derivation of this coverage model.

19.5 16.4 Level 3 — MSc Level: Probabilistic Models and Statistical Genomics

19.5.1 16.4.1 Hidden Markov Models — the Engine Inside HaplotypeCaller

A Pair-HMM models the alignment of a read to a candidate haplotype as a sequence of hidden states: Match (M), Insertion (I), Deletion (D). Each state has emission probabilities (how likely is this base given the state and quality score) and transition probabilities (gap-open/gap-extend penalties, derived from empirical indel rates).

This is literally the same HMM formalism used in speech recognition (phoneme alignment) and NLP (POS tagging) — applied to DNA. If a CS graduate already knows the Viterbi algorithm and the Forward-Backward algorithm from an NLP or ML course, HaplotypeCaller’s internals are immediately intelligible.

import numpy as np

def pair_hmm_forward(read, haplotype, match_prob=0.98, gap_open=0.02, gap_extend=0.4):
    """
    Simplified Pair-HMM forward algorithm — computes P(read | haplotype)
    by summing over ALL possible alignments (not just the best one, unlike Smith-Waterman).
    This full marginalization is what makes HaplotypeCaller's genotype likelihoods
    more robust than a single best-alignment approach.
    
    States: M (match/mismatch), I (insertion in read), D (deletion in read)
    Time: O(n*m*3) — three states per DP cell
    """
    n, m = len(read), len(haplotype)
    NEG_INF = -1e10
    
    # Log-space DP matrices for numerical stability (standard practice — avoids underflow
    # from multiplying many small probabilities)
    M = np.full((n+1, m+1), NEG_INF)
    I = np.full((n+1, m+1), NEG_INF)
    D = np.full((n+1, m+1), NEG_INF)
    
    M[0][0] = 0.0
    log_match = np.log(match_prob)
    log_mismatch = np.log((1-match_prob)/3)
    log_gap_open = np.log(gap_open)
    log_gap_extend = np.log(gap_extend)
    
    def logsumexp(*args):
        m_ = max(args)
        if m_ == NEG_INF: return NEG_INF
        return m_ + np.log(sum(np.exp(a - m_) for a in args))
    
    for i in range(1, n+1):
        for j in range(1, m+1):
            emit = log_match if read[i-1] == haplotype[j-1] else log_mismatch
            M[i][j] = emit + logsumexp(M[i-1][j-1], I[i-1][j-1], D[i-1][j-1])
            I[i][j] = logsumexp(M[i-1][j] + log_gap_open, I[i-1][j] + log_gap_extend)
            D[i][j] = logsumexp(M[i][j-1] + log_gap_open, D[i][j-1] + log_gap_extend)
    
    # Total probability = sum over all paths ending anywhere in the last row/column
    total_log_prob = logsumexp(M[n][m], I[n][m], D[n][m])
    return total_log_prob

read = "ACGTACGT"
hap_ref = "ACGTACGT"      # perfect match haplotype
hap_alt = "ACGTTCGT"      # 1 SNP haplotype (position 5: A->T)

p_ref = pair_hmm_forward(read, hap_ref)
p_alt = pair_hmm_forward(read, hap_alt)
print(f"log P(read | reference haplotype) = {p_ref:.4f}")
print(f"log P(read | alt haplotype)       = {p_alt:.4f}")
print(f"\nHaplotypeCaller sums these likelihoods across ALL reads covering a site")
print(f"to compute genotype likelihoods via Bayes' rule — this IS the pair-HMM step")
print(f"referenced in the GATK papers.")

Reference: Durbin R., Eddy S., Krogh A., Mitchison G. (1998). Biological Sequence Analysis: Probabilistic Models of Proteins and Nucleic Acids. Cambridge University Press — the canonical HMM-in-bioinformatics textbook.

19.5.2 16.4.2 Expectation-Maximization: Salmon/kallisto’s Secret Engine

When a read maps ambiguously to multiple transcript isoforms of the same gene, how do we assign fractional counts fairly? Expectation-Maximization (EM) — the same algorithm behind Gaussian Mixture Models in general ML.

import numpy as np

def em_transcript_quantification(read_transcript_compatibility, n_iterations=50):
    """
    Simplified EM for transcript abundance estimation (core of Salmon/kallisto/RSEM).
    
    read_transcript_compatibility: list of sets, each set = {transcript_ids this read is compatible with}
    
    E-step: given current abundance estimates, compute expected assignment of each 
            ambiguous read to each compatible transcript, weighted by relative abundance
    M-step: given expected assignments, re-estimate each transcript's abundance
    
    Converges to the Maximum Likelihood Estimate of transcript abundances.
    """
    transcripts = sorted(set.union(*read_transcript_compatibility))
    n_t = len(transcripts)
    t_idx = {t: i for i, t in enumerate(transcripts)}
    
    # Initialize: uniform abundance
    theta = np.ones(n_t) / n_t
    
    for iteration in range(n_iterations):
        # E-step: fractional assignment of each read to compatible transcripts
        expected_counts = np.zeros(n_t)
        for compat_set in read_transcript_compatibility:
            weights = np.array([theta[t_idx[t]] for t in compat_set])
            weights = weights / weights.sum()  # normalize among compatible transcripts
            for t, w in zip(compat_set, weights):
                expected_counts[t_idx[t]] += w
        
        # M-step: re-estimate abundance from expected counts
        theta = expected_counts / expected_counts.sum()
    
    return dict(zip(transcripts, theta))

# Example: gene with 3 isoforms; some reads map uniquely, some ambiguously (multi-mapping)
reads = [
    {'iso1'}, {'iso1'}, {'iso1'},           # 3 reads unique to isoform 1
    {'iso2'}, {'iso2'},                     # 2 reads unique to isoform 2
    {'iso3'},                               # 1 read unique to isoform 3
    {'iso1', 'iso2'}, {'iso1', 'iso2'},     # 2 ambiguous reads (iso1 or iso2)
    {'iso2', 'iso3'},                       # 1 ambiguous read (iso2 or iso3)
]

result = em_transcript_quantification(reads)
print("EM-estimated relative abundances (this is what becomes TPM in Salmon):")
for t, abundance in sorted(result.items()):
    print(f"  {t}: {abundance:.4f}")
print("\nNote how ambiguous reads are fractionally assigned proportional to")
print("each isoform's UNIQUE-read-supported abundance — exactly Salmon's approach.")

Reference: Dempster A.P., Laird N.M., Rubin D.B. (1977). “Maximum Likelihood from Incomplete Data via the EM Algorithm.” JRSS-B 39(1):1–38 — the foundational EM paper, applied to genomics by Li B., Dewey C.N. (2011) “RSEM” BMC Bioinformatics 12:323, and Patro R. et al. (2017) Salmon.

19.5.3 16.4.3 Generalized Linear Models: The Negative Binomial GLM Behind DESeq2

cat("
DESeq2's statistical model, in full:

For gene i, sample j:
   K_ij ~ NegativeBinomial(mean = mu_ij, dispersion = alpha_i)
   
   mu_ij = s_j * q_ij                     (size factor x true expression)
   log2(q_ij) = sum_r( x_jr * beta_ir )   (linear model on log2 scale — a GLM!)

Where:
   s_j    = per-sample size factor (library size normalization via median-of-ratios)
   x_jr   = design matrix (e.g., condition: treated=1/control=0)
   beta_ir = the log2 fold change we are estimating (the coefficient of interest)
   alpha_i = gene-specific dispersion, SHRUNK toward a fitted mean-dispersion trend
             across all genes (this shrinkage is empirical Bayes — borrowing strength
             across genes, exactly like variance shrinkage in limma/edgeR)

This is precisely a Generalized Linear Model with:
   - Negative Binomial family (variance = mu + alpha*mu^2, i.e. quadratic mean-variance)
   - log link function
   - Empirical Bayes shrinkage on the dispersion parameter (not on beta itself, 
     though lfcShrink() via apeglm ALSO shrinks beta with a heavy-tailed prior)
")
#> 
#> DESeq2's statistical model, in full:
#> 
#> For gene i, sample j:
#>    K_ij ~ NegativeBinomial(mean = mu_ij, dispersion = alpha_i)
#>    
#>    mu_ij = s_j * q_ij                     (size factor x true expression)
#>    log2(q_ij) = sum_r( x_jr * beta_ir )   (linear model on log2 scale — a GLM!)
#> 
#> Where:
#>    s_j    = per-sample size factor (library size normalization via median-of-ratios)
#>    x_jr   = design matrix (e.g., condition: treated=1/control=0)
#>    beta_ir = the log2 fold change we are estimating (the coefficient of interest)
#>    alpha_i = gene-specific dispersion, SHRUNK toward a fitted mean-dispersion trend
#>              across all genes (this shrinkage is empirical Bayes — borrowing strength
#>              across genes, exactly like variance shrinkage in limma/edgeR)
#> 
#> This is precisely a Generalized Linear Model with:
#>    - Negative Binomial family (variance = mu + alpha*mu^2, i.e. quadratic mean-variance)
#>    - log link function
#>    - Empirical Bayes shrinkage on the dispersion parameter (not on beta itself, 
#>      though lfcShrink() via apeglm ALSO shrinks beta with a heavy-tailed prior)
# Fit a simple NB GLM manually with MASS::glm.nb to show what DESeq2 does under the hood
if (requireNamespace("MASS", quietly=TRUE)) {
  library(MASS)
  set.seed(1)
  n <- 200
  condition <- rep(c(0,1), each=n/2)
  true_lfc <- 1.5
  mu <- 50 * 2^(condition * true_lfc)
  counts <- rnbinom(n, mu=mu, size=5)  # size = 1/dispersion
  
  df <- data.frame(counts=counts, condition=condition)
  fit <- glm.nb(counts ~ condition, data=df)
  
  cat("Manual NB GLM fit (single gene):\n")
  print(summary(fit)$coefficients)
  cat(sprintf("\nTrue log2FC used in simulation: %.2f\n", true_lfc))
  cat(sprintf("Estimated log2FC (beta/log(2)): %.2f\n", coef(fit)[2]/log(2)))
  cat("\nThis IS conceptually what DESeq2 does per-gene, at scale, with\n")
  cat("additional empirical-Bayes dispersion shrinkage across ~20,000 genes.\n")
} else {
  cat("Install MASS: install.packages('MASS')\n")
}
#> Manual NB GLM fit (single gene):
#>              Estimate Std. Error  z value     Pr(>|z|)
#> (Intercept) 3.9592882 0.04831977 81.93929 0.000000e+00
#> condition   0.9507161 0.06747265 14.09039 4.351403e-45
#> 
#> True log2FC used in simulation: 1.50
#> Estimated log2FC (beta/log(2)): 1.37
#> 
#> This IS conceptually what DESeq2 does per-gene, at scale, with
#> additional empirical-Bayes dispersion shrinkage across ~20,000 genes.

Reference: McCullagh P., Nelder J.A. (1989). Generalized Linear Models, 2nd ed. Chapman & Hall — the canonical GLM reference; Love M.I. et al. (2014). Genome Biology 15:550 (DESeq2’s specific NB-GLM + shrinkage formulation).

19.5.4 16.4.4 Graph Algorithms: De Bruijn Graphs for Assembly

De Bruijn graph assembly turns genome assembly into an Eulerian path problem — find a path visiting every edge exactly once. This is a classical CS graph theory problem (Eulerian vs. Hamiltonian paths), and its application to genome assembly (rather than the OLC/overlap-based Hamiltonian-path formulation used by Sanger-era assemblers) was the key algorithmic innovation that made short-read assembly computationally tractable.

Why Eulerian, not Hamiltonian? Finding a Hamiltonian path (visit every node once) is NP-hard. Finding an Eulerian path (visit every edge once) is solvable in linear time (Hierholzer’s algorithm). By representing each k-mer as an edge (not a node) in a de Bruijn graph, genome assembly transforms from an intractable to a tractable problem — this reformulation (Pevzner, Tang, Waterman, 2001) is one of the most important algorithmic insights in computational biology.

def build_de_bruijn_graph(reads, k=4):
    """
    Build a de Bruijn graph from sequencing reads.
    Nodes = (k-1)-mers; Edges = k-mers (each k-mer connects its prefix (k-1)-mer 
    to its suffix (k-1)-mer).
    Time: O(total sequence length)
    """
    graph = {}
    for read in reads:
        for i in range(len(read) - k + 1):
            kmer = read[i:i+k]
            prefix, suffix = kmer[:-1], kmer[1:]
            graph.setdefault(prefix, []).append(suffix)
    return graph

def find_eulerian_path(graph):
    """
    Hierholzer's algorithm for Eulerian path — O(E) time (E = number of edges).
    This IS how SPAdes/Velvet/ABySS conceptually assemble contigs from k-mers,
    though production tools add extensive error-correction and graph simplification.
    """
    graph = {k: list(v) for k, v in graph.items()}  # copy
    
    # Find start node (out_degree - in_degree = 1, or any node if graph is Eulerian circuit)
    out_deg = {n: len(v) for n, v in graph.items()}
    in_deg = {}
    for v in graph.values():
        for node in v:
            in_deg[node] = in_deg.get(node, 0) + 1
    
    start = next(iter(graph))
    for node in out_deg:
        if out_deg[node] - in_deg.get(node, 0) == 1:
            start = node
            break
    
    # Hierholzer's algorithm
    stack, path = [start], []
    while stack:
        node = stack[-1]
        if graph.get(node):
            stack.append(graph[node].pop())
        else:
            path.append(stack.pop())
    return path[::-1]

# Simulate reads from a genome with a repeat (classic assembly challenge)
genome = "ACGTACGTAGCTAGCTACGT"  # contains repeated "ACGT" — a real assembly complexity
reads = [genome[i:i+8] for i in range(0, len(genome)-7, 3)]  # simulate overlapping reads
print(f"Original genome: {genome}")
print(f"Simulated reads: {reads}")

dbg = build_de_bruijn_graph(reads, k=5)
path = find_eulerian_path(dbg)
assembled = path[0] + ''.join(node[-1] for node in path[1:])
print(f"\nAssembled sequence: {assembled}")
print(f"Note: repeats (like 'ACGT' appearing twice) can create AMBIGUOUS graph paths —")
print(f"this is precisely why repetitive regions are the hardest part of genome assembly,")
print(f"and why long reads (Nanopore/PacBio) that SPAN repeats are so valuable (Section 8).")

Reference: Pevzner P.A., Tang H., Waterman M.S. (2001). “An Eulerian path approach to DNA fragment assembly.” PNAS 98(17):9748–9753 — the paper that introduced de Bruijn graph assembly to genomics. [Link]

19.6 16.5 Level 4 — Pre-Doctoral: Advanced Statistical & ML Machinery

19.6.1 16.5.1 Variational Inference and Deep Generative Models: scVI

scVI (single-cell Variational Inference) models scRNA-seq counts with a Variational Autoencoder (VAE): an encoder network maps each cell’s expression vector to a latent Gaussian distribution (mean, variance), a decoder network maps latent samples back to a Zero-Inflated Negative Binomial (ZINB) distribution over observed counts, and the whole model is trained by maximizing the Evidence Lower Bound (ELBO) — exactly the same variational inference machinery from Kingma & Welling’s original VAE paper (2013), adapted to genomics’ zero-inflated, overdispersed count data. Batch effects are corrected by conditioning the decoder on a batch covariate, forcing the latent space to be batch-invariant.

The ELBO being optimized: \[\mathcal{L}(\theta,\phi;x) = \mathbb{E}_{q_\phi(z|x)}[\log p_\theta(x|z)] - D_{KL}(q_\phi(z|x) \| p(z))\]

  • First term: reconstruction likelihood (how well can we regenerate the observed counts from the latent code)
  • Second term: KL divergence regularizer (keeps the latent space close to a standard normal prior — this is what enables generative sampling and smooth interpolation)
# Conceptual PyTorch sketch of the scVI encoder/decoder architecture
# (For production, always use the scvi-tools package: pip install scvi-tools)
import torch
import torch.nn as nn

class EncoderVAE(nn.Module):
    """Maps gene expression -> latent distribution parameters (mu, log_var)."""
    def __init__(self, n_genes, n_latent=10, n_hidden=128):
        super().__init__()
        self.net = nn.Sequential(
            nn.Linear(n_genes, n_hidden), nn.ReLU(), nn.BatchNorm1d(n_hidden),
            nn.Linear(n_hidden, n_hidden), nn.ReLU(), nn.BatchNorm1d(n_hidden)
        )
        self.mu_layer  = nn.Linear(n_hidden, n_latent)
        self.var_layer = nn.Linear(n_hidden, n_latent)
    
    def forward(self, x):
        h = self.net(torch.log1p(x))  # log1p transform of raw counts
        mu = self.mu_layer(h)
        log_var = self.var_layer(h)
        return mu, log_var

class DecoderZINB(nn.Module):
    """Maps latent code (+ batch) -> Zero-Inflated Negative Binomial parameters."""
    def __init__(self, n_latent, n_genes, n_batch, n_hidden=128):
        super().__init__()
        self.net = nn.Sequential(
            nn.Linear(n_latent + n_batch, n_hidden), nn.ReLU(),
            nn.Linear(n_hidden, n_hidden), nn.ReLU()
        )
        self.mean_layer = nn.Linear(n_hidden, n_genes)   # NB mean (softmax normalized)
        self.dropout_layer = nn.Linear(n_hidden, n_genes)  # zero-inflation logits
    
    def forward(self, z, batch_onehot):
        h = self.net(torch.cat([z, batch_onehot], dim=1))
        px_scale = torch.softmax(self.mean_layer(h), dim=-1)  # relative expression
        px_dropout = self.dropout_layer(h)                     # dropout probability logits
        return px_scale, px_dropout

def reparameterize(mu, log_var):
    """Reparameterization trick — makes sampling differentiable for backprop."""
    std = torch.exp(0.5 * log_var)
    eps = torch.randn_like(std)
    return mu + eps * std

def elbo_loss(x, px_scale, px_dropout, library_size, mu, log_var, theta):
    """
    Simplified ELBO: ZINB reconstruction likelihood - KL divergence to N(0,I) prior.
    theta = learned per-gene dispersion parameter.
    """
    kl_div = -0.5 * torch.sum(1 + log_var - mu.pow(2) - log_var.exp(), dim=1)
    # (Full ZINB negative log-likelihood omitted for brevity — see scvi-tools source)
    return kl_div.mean()

print("scVI architecture sketch defined. In practice:")
print("  import scvi")
print("  scvi.model.SCVI.setup_anndata(adata, batch_key='batch')")
print("  model = scvi.model.SCVI(adata, n_latent=10)")
print("  model.train()")
print("  latent = model.get_latent_representation()  # batch-corrected embedding")

Reference: Lopez R. et al. (2018). “Deep generative modeling for single-cell transcriptomics.” Nature Methods 15:1053–1058 (scVI); Kingma D.P., Welling M. (2013). “Auto-Encoding Variational Bayes.” arXiv:1312.6114 (the foundational VAE paper).

19.6.2 16.5.2 Graph Neural Networks in Genomics

Modern applications extend GNNs (message-passing neural networks) to: (1) protein structure prediction — AlphaFold2’s Evoformer uses attention over a residue-pair graph; (2) gene regulatory network inference from scRNA-seq/scATAC-seq; (3) spatial transcriptomics — modeling cell-cell interactions as a graph where edges = physical proximity.

The message-passing framework, generally: \[h_v^{(k+1)} = \text{UPDATE}\left(h_v^{(k)}, \text{AGGREGATE}\left(\{h_u^{(k)} : u \in \mathcal{N}(v)\}\right)\right)\]

AlphaFold2’s Evoformer block is a specialized instance of this: nodes = residues, and the “graph” is implicitly the full pairwise residue-residue attention matrix, iteratively refined jointly with a Multiple Sequence Alignment (MSA) representation.

Reference: Jumper J. et al. (2021). “Highly accurate protein structure prediction with AlphaFold.” Nature 596:583–589. [Link]

# Minimal GNN message-passing layer sketch (PyTorch Geometric style, for spatial transcriptomics)
import torch
import torch.nn as nn

class SimpleGNNLayer(nn.Module):
    """
    One message-passing layer for a cell-cell interaction graph
    (nodes = cells with expression vectors, edges = spatial proximity in Visium/Xenium data).
    """
    def __init__(self, in_dim, out_dim):
        super().__init__()
        self.message_fn = nn.Linear(in_dim * 2, out_dim)
        self.update_fn  = nn.Linear(in_dim + out_dim, out_dim)
    
    def forward(self, node_features, edge_index):
        """
        node_features: [n_cells, in_dim]
        edge_index: [2, n_edges] — (source, target) pairs from spatial adjacency
        """
        src, dst = edge_index
        # Compute messages along each edge
        messages = self.message_fn(torch.cat([node_features[src], node_features[dst]], dim=1))
        # Aggregate messages per target node (sum aggregation)
        aggregated = torch.zeros(node_features.size(0), messages.size(1))
        aggregated.index_add_(0, dst, messages)
        # Update node representations
        updated = self.update_fn(torch.cat([node_features, aggregated], dim=1))
        return torch.relu(updated)

print("GNN layer defined for spatial/cell-cell-interaction graphs.")
print("Production use: squidpy.gr.spatial_neighbors() builds the graph,")
print("then a GNN (e.g., via PyTorch Geometric) can predict niche-level gene programs.")

19.6.3 16.5.3 Population Genetics Theory: Hardy-Weinberg, F_ST, and Coalescent Theory

# Hardy-Weinberg Equilibrium — the null model against which ALL population genetics deviates
# p^2 (AA) + 2pq (Aa) + q^2 (aa) = 1, where p+q=1

p <- seq(0, 1, 0.01)
q <- 1 - p
hwe_df <- data.frame(p=p, AA=p^2, Aa=2*p*q, aa=q^2)
hwe_long <- pivot_longer(hwe_df, cols=c(AA,Aa,aa), names_to="Genotype", values_to="Frequency")

ggplot(hwe_long, aes(x=p, y=Frequency, color=Genotype)) +
  geom_line(linewidth=1.3) +
  scale_color_manual(values=c("AA"="#e74c3c","Aa"="#f39c12","aa"="#3498db")) +
  labs(title="Hardy-Weinberg Equilibrium Genotype Frequencies",
       subtitle="Deviation from HWE (chi-square test) flags genotyping errors or population stratification",
       x="Allele frequency p (of allele A)", y="Genotype Frequency") +
  theme_minimal(base_size=12)

# Hardy-Weinberg equilibrium chi-square test — a standard GWAS QC step
from scipy import stats

def hwe_chisq_test(n_AA, n_Aa, n_aa):
    """
    Test whether observed genotype counts deviate from HWE.
    Used as a standard QC filter in GWAS: SNPs failing HWE (p < 1e-6 in controls)
    are typically excluded as likely genotyping artifacts.
    """
    n_total = n_AA + n_Aa + n_aa
    p = (2*n_AA + n_Aa) / (2*n_total)  # allele frequency of A
    q = 1 - p
    
    expected_AA = p**2 * n_total
    expected_Aa = 2*p*q * n_total
    expected_aa = q**2 * n_total
    
    observed = [n_AA, n_Aa, n_aa]
    expected = [expected_AA, expected_Aa, expected_aa]
    
    chi2, p_value = stats.chisquare(observed, expected, ddof=1)  # 1 dof (p estimated from data)
    return chi2, p_value

# Example: SNP genotype counts from a GWAS cohort
chi2, pval = hwe_chisq_test(n_AA=450, n_Aa=430, n_aa=120)
print(f"Chi-square statistic: {chi2:.3f}")
print(f"P-value: {pval:.4f}")
print(f"HWE QC: {'PASS' if pval > 1e-6 else 'FAIL — likely genotyping error, exclude SNP'}")

F_ST — quantifying population differentiation:

def fst_hudson(p1, p2, n1, n2):
    """
    Hudson's F_ST estimator — measures genetic differentiation between two populations.
    F_ST ranges 0 (no differentiation) to 1 (complete differentiation/fixation).
    
    p1, p2 = allele frequencies in population 1 and 2
    n1, n2 = sample sizes
    
    Used in: population structure analysis, selection scans (high F_ST = candidate
    for local adaptation), ancestry-informative marker selection.
    """
    numerator = (p1 - p2)**2 - (p1*(1-p1))/(n1-1) - (p2*(1-p2))/(n2-1)
    denominator = p1*(1-p2) + p2*(1-p1)
    return numerator / denominator if denominator != 0 else 0

# Example: a SNP with very different frequencies between two populations (candidate for selection)
fst = fst_hudson(p1=0.85, p2=0.15, n1=200, n2=200)
print(f"F_ST = {fst:.4f}")
print(f"Interpretation: {'High differentiation (candidate selection signal)' if fst > 0.15 else 'Low/moderate differentiation'}")

References: Hardy G.H. (1908). “Mendelian proportions in a mixed population.” Science 28(706):49–50; Weinberg W. (1908) [German, independently]; Hudson R.R., Slatkin M., Maddison W.P. (1992). “Estimation of levels of gene flow from DNA sequence data.” Genetics 132(2):583–589 (F_ST estimator).

19.6.4 16.5.4 Coalescent Theory — The Statistical Foundation of Population Genomics

Coalescent theory (Kingman, 1982) models the genealogy of a sample of genes backward in time, asking: how many generations until two randomly chosen lineages find their most recent common ancestor (MRCA)? For a population of size N, the expected time to coalescence for two lineages is N generations, and the full genealogy of a sample of n lineages follows the Kingman coalescent — a random binary tree with exponentially distributed branch lengths.

This is the theoretical backbone of: demographic inference (PSMC, MSMC — inferring historical population size from genome-wide heterozygosity patterns), ARG (Ancestral Recombination Graph) inference (tsinfer, Relate), and simulation tools (msprime, SLiM) used to generate null distributions for selection tests.

Reference: Kingman J.F.C. (1982). “The coalescent.” Stochastic Processes and their Applications 13(3):235–248; Li H., Durbin R. (2011). “Inference of human population history from individual whole-genome sequences.” Nature 475:493–496 (PSMC). [Link]

import numpy as np

def simulate_coalescent(n_lineages, effective_pop_size=10000):
    """
    Simulate coalescent times for a sample under the standard Kingman coalescent.
    Time to next coalescence (when k lineages remain) ~ Exponential(rate = C(k,2)/N)
    where C(k,2) = k*(k-1)/2 is the number of possible pairs that could coalesce.
    
    This IS the null model used to compute expected genetic diversity, 
    Tajima's D, and to calibrate demographic inference methods.
    """
    times = []
    k = n_lineages
    total_time = 0
    N = effective_pop_size
    
    while k > 1:
        rate = (k * (k-1) / 2) / N   # coalescence rate with k lineages
        t = np.random.exponential(1/rate)
        total_time += t
        times.append((k, t, total_time))
        k -= 1
    
    return times

np.random.seed(42)
coal_times = simulate_coalescent(n_lineages=10, effective_pop_size=10000)
print("Coalescent simulation (n=10 lineages, Ne=10,000):")
print(f"{'Lineages':>10} {'Wait time':>12} {'Cumulative':>12}")
for k, t, cum in coal_times:
    print(f"{k:>10} {t:>12.1f} {cum:>12.1f}")
print(f"\nTime to Most Recent Common Ancestor (MRCA): {coal_times[-1][2]:.0f} generations")
print("Note: coalescence accelerates as fewer lineages remain (more pairs -> higher rate")
print("early when there are lots of lineages, but the FIRST coalescent event actually")
print("happens fastest because C(k,2) is largest when k is largest).")

19.6.5 16.5.5 Causal Inference in Genomics: Mendelian Randomization

Mendelian Randomization (MR) exploits the random assortment of alleles at conception (Mendel’s laws) as a natural randomized controlled trial, using a genetic variant robustly associated with an exposure (e.g., LDL cholesterol) as an instrumental variable to estimate the variant’s causal effect on an outcome (e.g., coronary artery disease) — avoiding the confounding that plagues purely observational epidemiology.

Three core IV assumptions (borrowed directly from econometrics’ instrumental variables framework):
1. Relevance: the genetic variant is robustly associated with the exposure
2. Independence: the variant is not associated with confounders
3. Exclusion restriction: the variant affects the outcome only through the exposure, not via any other pathway

Reference: Davey Smith G., Ebrahim S. (2003). “‘Mendelian randomization’: can genetic epidemiology contribute to understanding environmental determinants of disease?” International Journal of Epidemiology 32(1):1–22.

# Two-Sample Mendelian Randomization — R workflow using the TwoSampleMR package
# install.packages("TwoSampleMR") or devtools::install_github("MRCIEU/TwoSampleMR")

# library(TwoSampleMR)
# 
# # Step 1: Extract instrument SNPs (genome-wide significant for exposure) from a GWAS
# exposure_dat <- extract_instruments(outcomes = 'ieu-a-2')  # e.g., LDL cholesterol GWAS
# 
# # Step 2: Extract the same SNPs' effects on the outcome GWAS
# outcome_dat <- extract_outcome_data(snps = exposure_dat$SNP, outcomes = 'ieu-a-7')  # CAD GWAS
# 
# # Step 3: Harmonize effect alleles between exposure and outcome datasets
# dat <- harmonise_data(exposure_dat, outcome_dat)
# 
# # Step 4: Run MR — Inverse Variance Weighted (IVW) is the primary method
# mr_results <- mr(dat, method_list = c("mr_ivw", "mr_egger_regression", "mr_weighted_median"))
# 
# # MR-Egger intercept test checks for horizontal pleiotropy (violation of exclusion restriction)
# pleiotropy_test <- mr_pleiotropy_test(dat)

cat("Mendelian Randomization workflow (see commented R code above for TwoSampleMR usage).\n")
cat("Key sensitivity analyses every rigorous MR study reports:\n")
cat("  1. IVW (primary estimate, assumes no pleiotropy)\n")
cat("  2. MR-Egger (allows pleiotropy, less statistical power)\n")
cat("  3. Weighted median (robust if <50% of instruments are invalid)\n")
cat("  4. MR-PRESSO (detects and corrects for outlier/pleiotropic SNPs)\n")

19.7 16.6 Comprehensive Pre-Doctoral Question Bank (Increasing Difficulty)

[L1] What is the time complexity of building a hash table index over all k-mers in a genome of length n, and what is the average lookup time?
Answer: Building: O(n) — one pass, each insertion O(1) amortized. Lookup: O(1) average case (assuming good hash distribution and O(k) string comparison to confirm exact match, so more precisely O(k) per lookup — but this is treated as O(1) since k is a small constant like 20-31).

[L2] Why is Smith-Waterman preferred over Needleman-Wunsch for read alignment against a full chromosome, and what is done in practice to make it computationally feasible at genome scale?
Answer: Smith-Waterman finds the best LOCAL alignment (ignoring non-matching flanks), which matches the biological reality that a 150bp read corresponds to a small window within a 250Mb chromosome — Needleman-Wunsch would force end-to-end alignment, which is meaningless here. However, full Smith-Waterman is O(n·m) — at genome scale (n=250M) this is computationally infeasible per read. In practice: BWA/Bowtie use FM-index seed-and-extend — find short exact seeds via O(m) backward search, then run Smith-Waterman only in the small local window around each seed candidate, not against the whole chromosome.

**[L3] Derive why RNA-seq counts show var(K) = mu + alpha*mu^2 under the Negative Binomial model, and explain the biological interpretation of each term.**
Answer: The NB distribution arises as a Poisson-Gamma mixture: if K|lambda ~ Poisson(lambda) and lambda ~ Gamma(shape=1/alpha, scale=alphamu), then marginally K ~ NB(mean=mu, dispersion=alpha) with Var(K) = mu + alphamu^2. Biologically: the “mu” term is pure Poisson sampling noise (shot noise from random RNA fragment capture — irreducible even with a “perfect” experiment). The “alpha*mu^2” term captures biological variability between replicate samples (individual-to-individual variation in true expression) — this is what a Poisson model, which assumes var=mean, completely misses. As mu grows, biological variance dominates technical (Poisson) variance, which is why highly expressed genes show more absolute variability across replicates even though their relative (CV) variability may decrease.

[L4] In a Pair-HMM used by HaplotypeCaller, why is the computation done in log-space rather than directly with probabilities, and what numerical algorithm is needed to correctly sum probabilities in log-space?
Answer: Multiplying many small probabilities (each read-base emission probability, e.g., ~0.001 for a mismatch) causes floating-point underflow after even a few dozen multiplications — a 150bp read’s naive joint probability could require representing numbers smaller than 10⁻³⁰⁰, beyond double-precision float range. Working in log-space converts products into sums (numerically stable), but this creates a new problem: how do you compute log(a+b) when you only have log(a) and log(b)? Naively exponentiating back would reintroduce underflow. The solution is the log-sum-exp trick: log(a+b) = log(a) + log(1 + exp(log(b)-log(a))), computed by factoring out the larger term first (max(log a, log b)) to keep the exponent argument close to zero, preserving numerical precision. This exact technique reappears throughout ML (softmax computation, Viterbi/Forward-Backward algorithms, variational inference ELBO computation).

[Pre-doctoral / Research] Design a statistical test to detect selection at a locus using coalescent theory as the null model, and describe what patterns in the data would constitute evidence for positive selection versus balancing selection.
Answer: Under neutral coalescent theory (Kingman coalescent, constant population size), the site frequency spectrum (SFS) — the distribution of how many samples carry the derived allele at each variable site — follows a predictable shape (proportional to 1/i for i copies of the derived allele, under the standard neutral model). Tajima’s D compares two estimators of the population mutation parameter theta: one based on the number of segregating sites (theta_W, sensitive to rare variants) and one based on average pairwise heterozygosity (theta_pi, sensitive to intermediate-frequency variants). Under neutrality, D≈0. Positive selection (a selective sweep) reduces diversity around the selected site and skews the SFS toward rare variants (excess of low-frequency variants as new mutations arise on the sweeping haplotype after the sweep) → Tajima’s D significantly negative. Balancing selection maintains multiple haplotypes at intermediate frequency longer than neutral expectation → excess of intermediate-frequency variants → Tajima’s D significantly positive. A rigorous test requires simulating the null distribution of D under a realistic demographic model (via msprime/SLiM coalescent simulation, not just the simplistic constant-N assumption) because non-selective demographic events (population bottlenecks, expansions) can also skew Tajima’s D and are frequently confounded with selection signals — this is why modern selection scans (e.g., iHS, XP-EHH combined with demographic-aware simulation) are preferred over Tajima’s D alone in publication-quality analyses.

19.8 16.7 Summary Table: Concept Bridge (CS Course → Bioinformatics Application)

bridge <- data.frame(
  CS_Course_Concept = c("Hash tables", "Dynamic programming (edit distance)",
                          "Suffix trees/arrays", "Graph theory (Eulerian paths)",
                          "Hidden Markov Models (NLP)", "EM algorithm (GMMs)",
                          "Generalized Linear Models", "Variational Autoencoders",
                          "Graph Neural Networks", "A* / Dijkstra shortest path",
                          "Bloom filters", "MapReduce/distributed computing",
                          "Regular expressions/automata", "Bayesian inference",
                          "Compression algorithms (Huffman/LZ)"),
  Bioinformatics_Application = c(
    "K-mer indexing (Jellyfish, minimap2 seeding)",
    "Sequence alignment (Needleman-Wunsch, Smith-Waterman, BWA-MEM extension)",
    "Genome indexing for O(m) search (BWA's FM-index)",
    "De Bruijn graph genome assembly (SPAdes, Velvet)",
    "Variant calling (GATK Pair-HMM), gene prediction (exon/intron state)",
    "Transcript quantification (Salmon, kallisto, RSEM)",
    "Differential expression (DESeq2, edgeR negative binomial GLM)",
    "Batch-corrected embeddings for scRNA-seq (scVI, totalVI)",
    "Protein structure prediction (AlphaFold2), spatial cell-cell interactions",
    "Not directly common, but used in some haplotype phasing algorithms",
    "Fast approximate membership testing in metagenomic classifiers (Kraken2, Bloom-filter-based)",
    "Distributed variant calling (GATK Spark, Hail for population genomics)",
    "Sequence motif matching, restriction site finding, regex-based annotation parsing",
    "ACMG variant classification, GWAS fine-mapping, phylogenetics",
    "BAM/CRAM compression (reference-based delta encoding in CRAM)"
  ),
  stringsAsFactors=FALSE
)
kable(bridge, caption="Bridging standard CS curriculum to bioinformatics — a map for CS graduates") %>%
  kable_styling(bootstrap_options=c("striped","hover","condensed"), full_width=TRUE, font_size=11) %>%
  column_spec(1, bold=TRUE, color="white", background="#2c3e50")
Bridging standard CS curriculum to bioinformatics — a map for CS graduates
CS_Course_Concept Bioinformatics_Application
Hash tables K-mer indexing (Jellyfish, minimap2 seeding)
Dynamic programming (edit distance) Sequence alignment (Needleman-Wunsch, Smith-Waterman, BWA-MEM extension)
Suffix trees/arrays Genome indexing for O(m) search (BWA’s FM-index)
Graph theory (Eulerian paths) De Bruijn graph genome assembly (SPAdes, Velvet)
Hidden Markov Models (NLP) Variant calling (GATK Pair-HMM), gene prediction (exon/intron state)
EM algorithm (GMMs) Transcript quantification (Salmon, kallisto, RSEM)
Generalized Linear Models Differential expression (DESeq2, edgeR negative binomial GLM)
Variational Autoencoders Batch-corrected embeddings for scRNA-seq (scVI, totalVI)
Graph Neural Networks Protein structure prediction (AlphaFold2), spatial cell-cell interactions
A* / Dijkstra shortest path Not directly common, but used in some haplotype phasing algorithms
Bloom filters Fast approximate membership testing in metagenomic classifiers (Kraken2, Bloom-filter-based)
MapReduce/distributed computing Distributed variant calling (GATK Spark, Hail for population genomics)
Regular expressions/automata Sequence motif matching, restriction site finding, regex-based annotation parsing
Bayesian inference ACMG variant classification, GWAS fine-mapping, phylogenetics
Compression algorithms (Huffman/LZ) BAM/CRAM compression (reference-based delta encoding in CRAM)

20 Complete Reference List

20.1 Primary Tool Papers

  • Andrews S. (2010). FastQC. Babraham Bioinformatics. [Link]
  • Cock P.J.A. et al. (2010). “The Sanger FASTQ file format.” Nucleic Acids Research 38(6):1767–1771. [Link]
  • Ewing B. & Green P. (1998). “Base-calling of automated sequencer traces using phred II.” Genome Research 8:186–194. [PubMed]
  • Chen S. et al. (2018). “fastp: an ultra-fast all-in-one FASTQ preprocessor.” Bioinformatics 34(17):i884–i890. [Link]
  • Ewels P. et al. (2016). “MultiQC: summarize analysis results for multiple tools and samples.” Bioinformatics 32(19):3047–3048. [Link]
  • Li H. (2013). “Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM.” arXiv:1303.3997. [arXiv]
  • Vasimuddin Md. et al. (2019). “Efficient Architecture-Aware Acceleration of BWA-MEM for Multicore Systems.” IPDPS. [arXiv]
  • Li H. (2018). “Minimap2: pairwise alignment for nucleotide sequences.” Bioinformatics 34(18):3094–3100. [Link]
  • Li H. et al. (2009). “The Sequence Alignment/Map format and SAMtools.” Bioinformatics 25(16):2078–2079. [Link]
  • DePristo M.A. et al. (2011). “A framework for variation discovery and genotyping using next-generation DNA sequencing data.” Nature Genetics 43:491–498. [Link]
  • Van der Auwera G.A. & O’Connor B.D. (2020). Genomics in the Cloud. O’Reilly Media.
  • Poplin R. et al. (2018). “Scaling accurate genetic variant discovery to tens of thousands of samples.” bioRxiv 201178. [bioRxiv]
  • Poplin R. et al. (2018). “A universal SNP and small-indel variant caller using deep neural networks.” Nature Biotechnology 36:983–987. [Link]
  • Benjamin D. et al. (2019). “Calling Somatic SNVs and Indels with Mutect2.” bioRxiv 861054. [bioRxiv]
  • Cibulskis K. et al. (2013). “Sensitive detection of somatic point mutations in impure and heterogeneous cancer samples.” Nature Biotechnology 31:213–219.
  • Chen X. et al. (2016). “Manta: rapid detection of structural variants and indels.” Bioinformatics 32(8):1220–1222.
  • Talevich E. et al. (2016). “CNVkit: Genome-Wide Copy Number Detection.” PLOS Computational Biology 12(4):e1004873. [Link]
  • Danecek P. et al. (2011). “The variant call format and VCFtools.” Bioinformatics 27(15):2156–2158.
  • McLaren W. et al. (2016). “The Ensembl Variant Effect Predictor.” Genome Biology 17:122. [Link]
  • Chen S. et al. (2024). “A genomic mutational constraint map using variation in 76,156 human genomes.” Nature 625:92–100 (gnomAD v4). [Link]
  • Landrum M.J. et al. (2018). “ClinVar: improving access to variant interpretations.” Nucleic Acids Research 46(D1):D1062–D1067.
  • Richards S. et al. (2015). “Standards and guidelines for the interpretation of sequence variants.” Genetics in Medicine 17:405–424. [Link]
  • Dobin A. et al. (2013). “STAR: ultrafast universal RNA-seq aligner.” Bioinformatics 29(1):15–21.
  • Patro R. et al. (2017). “Salmon provides fast and bias-aware quantification.” Nature Methods 14:417–419. [Link]
  • Love M.I., Huber W., Anders S. (2014). “Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2.” Genome Biology 15:550. [Link]
  • Robinson M.D., McCarthy D.J., Smyth G.K. (2010). “edgeR: a Bioconductor package for differential expression analysis.” Bioinformatics 26(1):139–140.
  • Subramanian A. et al. (2005). “Gene set enrichment analysis.” PNAS 102(43):15545–15550. [Link]
  • Liao Y., Smyth G.K., Shi W. (2014). “featureCounts.” Bioinformatics 30(7):923–930.
  • Zheng G.X.Y. et al. (2017). “Massively parallel digital transcriptional profiling of single cells.” Nature Communications 8:14049.
  • Kaminow B., Yunusov D., Dobin A. (2021). “STARsolo.” bioRxiv 2021.05.05.442755.
  • Wolf F.A., Angerer P., Theis F.J. (2018). “SCANPY: large-scale single-cell gene expression data analysis.” Genome Biology 19:15.
  • Luecken M.D. & Theis F.J. (2019). “Current best practices in single-cell RNA-seq analysis.” Molecular Systems Biology 15:e8746.
  • Hafemeister C., Satija R. (2019). “Normalization and variance stabilization of single-cell RNA-seq data using SCTransform.” Genome Biology 20:296.
  • McInnes L., Healy J., Melville J. (2018). “UMAP: Uniform Manifold Approximation and Projection.” arXiv:1802.03426.
  • Traag V.A., Waltman L., van Eck N.J. (2019). “From Louvain to Leiden.” Scientific Reports 9:5233.
  • Hao Y. et al. (2021). “Integrated analysis of multimodal single-cell data.” Cell 184(13):3573–3587.
  • Korsunsky I. et al. (2019). “Fast, sensitive and accurate integration of single-cell data with Harmony.” Nature Methods 16:1289–1296.
  • Bergen V. et al. (2020). “Generalizing RNA velocity to transient cell states through dynamical modeling.” Nature Biotechnology 38:1408–1414.
  • Jin S. et al. (2021). “Inference and analysis of cell-cell communication using CellChat.” Nature Communications 12:1088.
  • Efremova M. et al. (2020). “CellPhoneDB.” Nature Protocols 15:1484–1506.
  • Tirosh I. et al. (2016). “Dissecting the multicellular ecosystem of metastatic melanoma by single-cell RNA-seq.” Science 352:189–196.
  • Gao R. et al. (2021). “Delineating copy number and clonal substructure from single-cell transcriptomes.” Nature Biotechnology 39:599–608 (CopyKAT).
  • Büttner M. et al. (2021). “scCODA is a Bayesian model for compositional single-cell data analysis.” Nature Communications 12:6876.
  • Wherry E.J., Kurachi M. (2015). “Molecular and cellular insights into T cell exhaustion.” Nature Reviews Immunology 15:486–499.
  • Miller B.C. et al. (2019). “Subsets of exhausted CD8+ T cells differentially mediate tumor control.” Nature Immunology 20:326–336.
  • Di Tommaso P. et al. (2017). “Nextflow enables reproducible computational workflows.” Nature Biotechnology 35:316–319.
  • Ewels P.A. et al. (2020). “The nf-core framework for community-curated bioinformatics pipelines.” Nature Biotechnology 38:276–278.
  • Nurk S. et al. (2022). “The complete sequence of a human genome.” Science 376:44–53. [Link]
  • Wang Y. et al. (2021). “Nanopore sequencing technology, bioinformatics and applications.” Nature Biotechnology 39:1348–1365.
  • Wenger A.M. et al. (2019). “Accurate circular consensus long-read sequencing improves variant detection.” Nature Biotechnology 37:1155–1162.
  • Benjamini Y., Hochberg Y. (1995). “Controlling the False Discovery Rate.” JRSS-B 57(1):289–300. [Link]
  • Li M.M. et al. (2017). “Standards and Guidelines for the Interpretation and Reporting of Sequence Variants in Cancer.” JMD 19(1):4–23.
  • Burrows M., Wheeler D.J. (1994). “A block-sorting lossless data compression algorithm.” DEC Technical Report 124. [PDF]
  • Sanger F., Nicklen S., Coulson A.R. (1977). “DNA sequencing with chain-terminating inhibitors.” PNAS 74(12):5463–5467. [Link]
  • Lander E.S. et al. (2001). “Initial sequencing and analysis of the human genome.” Nature 409:860–921.
  • Wolock S.L., Lopez R., Klein A.M. (2019). “Scrublet: Computational Identification of Cell Doublets.” Cell Systems 8(4):281–291.
  • Aran D. et al. (2019). “Reference-based analysis of lung single-cell sequencing.” Nature Immunology 20:163–172 (SingleR).
  • Stuart T. et al. (2019). “Comprehensive Integration of Single-Cell Data.” Cell 177(7):1888–1902.
  • Lun A.T.L. et al. (2019). “EmptyDrops.” Genome Biology 20:63.
  • Cheng J. et al. (2023). “Accurate proteome-wide missense variant effect prediction with AlphaMissense.” Science 381:eadg7492.
  • Kircher M. et al. (2014). “A general framework for estimating the relative pathogenicity of human genetic variants.” Nature Genetics 46:310–315 (CADD).
  • Jaganathan K. et al. (2019). “Predicting Splicing from Primary Sequence with Deep Learning.” Cell 176:535–548 (SpliceAI).
  • Li M.M. et al. (2017). “Standards and Guidelines for Somatic Variant Interpretation.” JMD 19(1):4–23.

20.2 Part 16 — CS Foundations & Pre-Doctoral References

  • Watson J.D., Crick F.H.C. (1953). “Molecular structure of nucleic acids.” Nature 171:737–738. [Link]
  • Crick F. (1958). “On Protein Synthesis.” Symposia of the Society for Experimental Biology 12:138–163 (origin of “Central Dogma”).
  • Needleman S.B., Wunsch C.D. (1970). “A general method applicable to the search for similarities in the amino acid sequence of two proteins.” Journal of Molecular Biology 48(3):443–453.
  • Smith T.F., Waterman M.S. (1981). “Identification of common molecular subsequences.” Journal of Molecular Biology 147(1):195–197.
  • Manber U., Myers G. (1993). “Suffix arrays: a new method for on-line string searches.” SIAM Journal on Computing 22(5):935–948.
  • Nong G., Zhang S., Chen W.H. (2009). “Linear Suffix Array Construction by Almost Pure Induced-Sorting.” Data Compression Conference (DCC).
  • Ferragina P., Manzini G. (2000). “Opportunistic data structures with applications.” FOCS 2000 (the original FM-index paper).
  • Lander E.S., Waterman M.S. (1988). “Genomic mapping by fingerprinting random clones: a mathematical analysis.” Genomics 2(3):231–239.
  • Marçais G., Kingsford C. (2011). “A fast, lock-free approach for efficient parallel counting of occurrences of k-mers.” Bioinformatics 27(6):764–770. [Link]
  • Durbin R., Eddy S., Krogh A., Mitchison G. (1998). Biological Sequence Analysis: Probabilistic Models of Proteins and Nucleic Acids. Cambridge University Press.
  • Dempster A.P., Laird N.M., Rubin D.B. (1977). “Maximum Likelihood from Incomplete Data via the EM Algorithm.” JRSS-B 39(1):1–38.
  • Li B., Dewey C.N. (2011). “RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome.” BMC Bioinformatics 12:323.
  • McCullagh P., Nelder J.A. (1989). Generalized Linear Models, 2nd ed. Chapman & Hall.
  • Pevzner P.A., Tang H., Waterman M.S. (2001). “An Eulerian path approach to DNA fragment assembly.” PNAS 98(17):9748–9753. [Link]
  • Lopez R. et al. (2018). “Deep generative modeling for single-cell transcriptomics.” Nature Methods 15:1053–1058.
  • Kingma D.P., Welling M. (2013). “Auto-Encoding Variational Bayes.” arXiv:1312.6114. [arXiv]
  • Jumper J. et al. (2021). “Highly accurate protein structure prediction with AlphaFold.” Nature 596:583–589. [Link]
  • Hardy G.H. (1908). “Mendelian proportions in a mixed population.” Science 28(706):49–50.
  • Hudson R.R., Slatkin M., Maddison W.P. (1992). “Estimation of levels of gene flow from DNA sequence data.” Genetics 132(2):583–589.
  • Kingman J.F.C. (1982). “The coalescent.” Stochastic Processes and their Applications 13(3):235–248.
  • Li H., Durbin R. (2011). “Inference of human population history from individual whole-genome sequences.” Nature 475:493–496 (PSMC). [Link]
  • Davey Smith G., Ebrahim S. (2003). “‘Mendelian randomization’: can genetic epidemiology contribute to understanding environmental determinants of disease?” International Journal of Epidemiology 32(1):1–22.
  • Tajima F. (1989). “Statistical method for testing the neutral mutation hypothesis by DNA polymorphism.” Genetics 123(3):585–595 (Tajima’s D).
  • Kelleher J., Etheridge A.M., McVean G. (2016). “Efficient Coalescent Simulation and Genealogical Analysis for Large Sample Sizes.” PLOS Computational Biology 12(5):e1004842 (msprime).

Document compiled for bioinformatics and genomics graduate training. Commands and tool versions change — always verify against current official documentation before production use.

Version: 2024. Last updated: August 2026.

sessionInfo()
#> R version 4.6.1 (2026-06-24 ucrt)
#> Platform: x86_64-w64-mingw32/x64
#> Running under: Windows 11 x64 (build 26200)
#> 
#> Matrix products: default
#>   LAPACK version 3.12.1
#> 
#> locale:
#> [1] LC_COLLATE=English_India.utf8  LC_CTYPE=English_India.utf8   
#> [3] LC_MONETARY=English_India.utf8 LC_NUMERIC=C                  
#> [5] LC_TIME=English_India.utf8    
#> 
#> time zone: Asia/Calcutta
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats4    stats     graphics  grDevices utils     datasets  methods  
#> [8] base     
#> 
#> other attached packages:
#>  [1] MASS_7.3-65                 Biostrings_2.80.1          
#>  [3] XVector_0.52.0              Seurat_5.5.1               
#>  [5] SeuratObject_5.4.0          sp_2.2-3                   
#>  [7] org.Hs.eg.db_3.23.1         AnnotationDbi_1.74.0       
#>  [9] clusterProfiler_4.20.0      DESeq2_1.52.0              
#> [11] SummarizedExperiment_1.42.0 Biobase_2.72.0             
#> [13] MatrixGenerics_1.24.0       matrixStats_1.5.0          
#> [15] GenomicRanges_1.64.0        Seqinfo_1.2.0              
#> [17] IRanges_2.46.0              S4Vectors_0.50.1           
#> [19] BiocGenerics_0.58.1         generics_0.1.4             
#> [21] apeglm_1.34.0               ggrepel_0.9.8              
#> [23] gridExtra_2.3.1             RColorBrewer_1.1-3         
#> [25] kableExtra_1.4.1            knitr_1.51                 
#> [27] scales_1.4.0                tidyr_1.3.2                
#> [29] dplyr_1.2.1                 ggplot2_4.0.3              
#> 
#> loaded via a namespace (and not attached):
#>   [1] fs_2.1.0                spatstat.sparse_3.2-0   enrichplot_1.32.0      
#>   [4] httr_1.4.8              numDeriv_2016.8-1.1     tools_4.6.1            
#>   [7] sctransform_0.4.3       R6_2.6.1                uwot_0.2.4             
#>  [10] lazyeval_0.2.3          mgcv_1.9-4              withr_3.0.3            
#>  [13] progressr_1.0.0         cli_3.6.6               textshaping_1.0.5      
#>  [16] spatstat.explore_3.8-2  fastDummies_1.7.6       scatterpie_0.2.6       
#>  [19] labeling_0.4.3          sass_0.4.10             mvtnorm_1.4-2          
#>  [22] S7_0.2.2                spatstat.data_3.1-9     ggridges_0.5.7         
#>  [25] pbapply_1.7-4           systemfonts_1.3.2       yulab.utils_0.2.4      
#>  [28] gson_0.2.1              DOSE_4.6.0              svglite_2.2.2          
#>  [31] parallelly_1.48.0       bbmle_1.0.25.1          rstudioapi_0.19.0      
#>  [34] RSQLite_3.53.3          gridGraphics_0.5-1      ica_1.0-3              
#>  [37] spatstat.random_3.5-1   GO.db_3.23.1            Matrix_1.7-5           
#>  [40] abind_1.4-8             lifecycle_1.0.5         yaml_2.3.12            
#>  [43] qvalue_2.44.0           SparseArray_1.12.2      Rtsne_0.17             
#>  [46] grid_4.6.1              blob_1.3.0              promises_1.5.0         
#>  [49] crayon_1.5.3            bdsmatrix_1.3-7         miniUI_0.1.2           
#>  [52] ggtangle_0.1.2          lattice_0.22-9          cowplot_1.2.0          
#>  [55] KEGGREST_1.52.2         pillar_1.11.1           future.apply_1.20.2    
#>  [58] codetools_0.2-20        glue_1.8.1              ggiraph_0.9.6          
#>  [61] ggfun_0.2.1             spatstat.univar_3.2-0   fontLiberation_0.1.0   
#>  [64] data.table_1.18.4       vctrs_0.7.3             png_0.1-9              
#>  [67] treeio_1.36.1           spam_2.11-4             gtable_0.3.6           
#>  [70] emdbook_1.3.14          cachem_1.1.0            xfun_0.60              
#>  [73] S4Arrays_1.12.0         mime_0.13               coda_0.19-4.1          
#>  [76] survival_3.8-6          aisdk_1.4.12            fitdistrplus_1.2-6     
#>  [79] ROCR_1.0-12             nlme_3.1-169            ggtree_4.2.0           
#>  [82] bit64_4.8.4             fontquiver_0.2.1        RcppAnnoy_0.0.23       
#>  [85] bslib_0.12.0            irlba_2.3.7             KernSmooth_2.23-26     
#>  [88] otel_0.2.0              DBI_1.3.0               tidyselect_1.2.1       
#>  [91] processx_3.9.0          bit_4.6.0               compiler_4.6.1         
#>  [94] httr2_1.3.0             xml2_1.6.0              fontBitstreamVera_0.1.1
#>  [97] DelayedArray_0.38.2     plotly_4.12.1           lmtest_0.9-40          
#> [100] callr_3.8.0             rappdirs_0.3.4          goftest_1.2-3          
#> [103] stringr_1.6.0           digest_0.6.39           spatstat.utils_3.2-4   
#> [106] rmarkdown_2.31          htmltools_0.5.9         pkgconfig_2.0.3        
#> [109] fastmap_1.2.0           rlang_1.3.0             htmlwidgets_1.6.4      
#> [112] shiny_1.14.0            farver_2.1.2            jquerylib_0.1.4        
#> [115] zoo_1.9-0               jsonlite_2.0.0          BiocParallel_1.46.0    
#> [118] GOSemSim_2.38.3         magrittr_2.0.5          ggplotify_0.1.3        
#> [121] dotCall64_1.2           patchwork_1.3.2         Rcpp_1.1.2             
#> [124] ape_5.8-1               ggnewscale_0.5.2        gdtools_0.5.1          
#> [127] reticulate_1.46.0       stringi_1.8.9           plyr_1.8.9             
#> [130] parallel_4.6.1          listenv_1.0.0           deldir_2.0-4           
#> [133] splines_4.6.1           tensor_1.5.1            locfit_1.5-9.12        
#> [136] igraph_2.3.3            spatstat.geom_3.8-2     enrichit_0.2.1         
#> [139] RcppHNSW_0.7.0          reshape2_1.4.5          evaluate_1.0.5         
#> [142] tweenr_2.0.3            httpuv_1.6.17           RANN_2.6.2             
#> [145] purrr_1.2.2             polyclip_1.10-7         future_1.75.0          
#> [148] scattermore_1.2         ggforce_0.5.0           xtable_1.8-8           
#> [151] RSpectra_0.16-2         tidytree_0.4.8          tidydr_0.0.6           
#> [154] later_1.4.8             viridisLite_0.4.3       tibble_3.3.1           
#> [157] aplot_0.3.1             memoise_2.0.1           cluster_2.1.8.2        
#> [160] globals_0.19.1