What this is. A self-contained course that takes you from "what is a FASTQ file" to building and honestly evaluating deep-learning and generative models on genomic, single-cell, spatial, and histopathology data, and applying them to disease biology and drug discovery. Sixteen modules, roughly 200 000 words, every concept explained from first principles, with runnable code, comparison tables, worked exercises, and ten graded capstone projects. Nothing here says "see the literature" in place of an explanation.

1. Who this is for

The course is written for three readers at once, and each one needs a different route through it.

There is no prerequisite beyond high-school biology and the willingness to type commands into a terminal. Everything else is built in this course.

2. The shape of the course

Seven parts, sixteen modules. Each module is independently readable but written to be read in order.

Plus two reference files you will use constantly: GLOSSARY.md (400+ terms) and CHEATSHEETS.md (20 quick-reference sheets).

The dependency structure is not strictly linear — this map shows what each module actually needs before it makes sense:

3. How to use it

Read actively, with a terminal open. Every code block in this course is written to be run after changing paths. Reading code is not learning code. The single highest-yield habit: after each section, reproduce its main command or plot on data you downloaded yourself.

Do the exercises. Each module ends with 5–8 graded exercises (Warm-up / Core / Stretch) and real solutions. The Stretch exercises are where competence actually forms.

Work the capstones. Module 16 contains ten projects with named public datasets, deliverables, traps, and grading rubrics. A portfolio of three completed capstones is worth more in a job interview than any certificate.

Suggested pace. Module 16 contains a 12-month week-by-week plan in three intensities (10 h/week, 20 h/week, full-time). Rough guide:

Three study paths, if you want a shorter route to a specific job:

• Omics analyst: 1, 2, 3, 4, 5, 6, 9, 15 + capstones 1, 2, 3.
• AI-for-biology engineer: 2, 3, 5, 6, 9, 10, 11, 12, 15 + capstones 7, 8, 10.
• Translational / drug discovery: 2, 5, 9, 13, 14, 15 + capstones 6, 9, plus Module 8 if tissue imaging is in scope.

4. What this course assumes about honesty

Three convictions run through every module, because they are what separates results that hold from results that evaporate.

1. A model's accuracy is a claim about a population, not a dataset. Nearly every inflated result in computational biology comes from a split that leaked — patient, site, batch, scanner, relative, or time. Module 9 gives the full taxonomy of leakage; every later module applies it.
2. Biological variation lives between donors, not between cells. You can sequence a million cells from three mice and still have n = 3. Modules 5, 6, and 7 are explicit about which unit of replication each test is entitled to use.
3. A method's output is not evidence that the method was appropriate. Clustering always returns clusters, enrichment always returns pathways, docking always returns poses, and a generative model always returns molecules. Each module names the failure mode alongside the method.

5. Conventions used throughout

• inline code is a command, filename, function, or identifier you type literally.
• Code blocks are tagged by language: bash, python, r. Paths in angle brackets (<your_sample>) are placeholders you replace.
• Formulas are written in LaTeX and always followed by a sentence defining every symbol in words.
• Cross-references are by module name, e.g. "Module 9 (Machine Learning)".
• Common pitfalls sections are not filler. They are the accumulated cost of other people's wasted months.
• Tool recommendations state the version era they were written against. Tooling in this field turns over every 18–24 months; the reasoning for choosing a tool is the durable part, and that is what each comparison table is designed to teach.

6. Before you start: the 30-minute setup

Do this now, not later. Module 1 covers it in detail, but the minimum is:

# 1. Install a conda distribution (miniforge is the least painful choice)
#    macOS/Linux:
curl -L -O "https://github.com/conda-forge/miniforge/releases/latest/download/Miniforge3-$(uname)-$(uname -m).sh"
bash Miniforge3-$(uname)-$(uname -m).sh -b -p "$HOME/miniforge3"
"$HOME/miniforge3/bin/conda" init "$(basename "$SHELL")"
# restart your shell

# 2. Make channel order permanent and correct (this ordering matters — see Module 1)
conda config --add channels defaults
conda config --add channels bioconda
conda config --add channels conda-forge
conda config --set channel_priority strict

# 3. A first environment to confirm everything works
conda create -n bioai-core -y python=3.11 numpy pandas scipy matplotlib seaborn \
    scikit-learn jupyterlab biopython pysam samtools
conda activate bioai-core
python -c "import pandas, Bio, pysam, sklearn; print('ok')"
samtools --version | head -1

If those last two lines print without error, you are ready for Module 1. If you are on Windows, install WSL2 first and work inside the Linux filesystem — not /mnt/c — or you will spend the course fighting file permissions and line endings.

7. A note on compute

Most of this course runs on a laptop. The exceptions, and the honest minimum for each:

Free options that genuinely work for learning: Google Colab (free T4 GPU, enough for Modules 9–12 exercises), Galaxy (no-code NGS pipelines), Terra and the UCSC/Ensembl public resources. Module 1 and Module 15 cover HPC and cloud properly, including how not to spend money by accident.

Start with Module 1 — The Working Environment and Toolkit.

Module 1 — The Working Environment and Toolkit

In one paragraph. This module sets up the physical and software environment you will use for every later module: where to actually run computational biology work (laptop, shared HPC cluster, or cloud), the minimum set of command-line skills that let you move, inspect, and filter large text files without opening them, and a reproducible way to install scientific software using conda/mamba so that a pipeline you build today still runs in a year. By the end you will have four ready-to-use environment specifications covering genomics, single-cell omics, digital pathology, and deep learning, and you will understand why each tool lives where it does in those files.

Prerequisites: basic comfort using a computer (any OS); no prior programming or Linux experience assumed. You will be able to: - Decide, for a given analysis, whether a laptop, an HPC cluster, or a cloud instance is the right place to run it - Navigate a Unix filesystem and chain commands with pipes to answer real questions about sequencing files - Use grep, awk, sed, sort, cut, find, and xargs together to filter and reshape text-based bioinformatics files - Keep a long job running after disconnecting from a remote machine, using tmux/screen or nohup - Move files reliably between your laptop and a remote server with scp and rsync - Read and change Unix file permissions and understand why shared HPC systems enforce them - Build and install conda/mamba environments with correct channel priority, and explain the pip-inside-conda rule - Stand up four complete, working environment files for genomics, single-cell, pathology, and deep-learning work

Time: 4-6 hours (longer if this is your first exposure to a Unix shell)

1.1 How computational biology work is actually organised

Every computational biology project eventually needs compute (a machine to run code), storage (a place to keep files), and a way to get results out. The question "where do I run this?" has three real answers, and picking the wrong one wastes more time than any software bug.

Laptop. Your own machine, with a few to tens of gigabytes of RAM and a handful of CPU cores. Good for: writing and debugging code, working on small or subsampled datasets, exploratory plotting, anything involving a GUI (viewing a BAM file in a genome browser, reading an image in a slide viewer). Bad for: aligning hundreds of whole-genome FASTQ files, training a neural network on millions of cells, anything that needs more memory than your laptop has — a single human genome alignment index for bwa alone needs roughly 4-6 GB of RAM just to load, before you add the read data.

HPC (high-performance computing) cluster. A shared pool of machines ("nodes") managed by your institution, accessed by logging in to a "login node" and submitting jobs to a scheduler (commonly Slurm, sometimes PBS or LSF) that runs your job on a "compute node" when resources are free. Good for: large batch jobs (aligning hundreds of samples, variant calling across a cohort), anything that needs many CPU cores or large RAM (some single-cell analyses need 256 GB+ for big atlases), long-running jobs (multi-day phylogenetics or assembly runs), work that must stay inside an institutional or hospital network for data governance reasons (patient genomes rarely leave a firewall). Costs you nothing extra per hour (usually) but queueing for a busy cluster can mean waiting hours or days, and software you want may not be installed system-wide — which is exactly why you bring your own conda environment.

Cloud (AWS, Google Cloud, Azure, or specialised platforms such as DNAnexus or Terra). You rent exactly the machine you need, for exactly as long as you need it, and pay per hour. Good for: elastic burst capacity (you need 500 CPU-hours this week and zero next week), GPU access for deep learning when your institution has none or they are all booked, reproducible pipelines shared with external collaborators, production services. The real cost is not just the hourly rate — it is data egress fees (moving data out of cloud storage can cost more than the compute did) and the discipline required to shut instances down. A p3.2xlarge GPU instance left running idle over a weekend by mistake is a genuinely common way labs burn through grant money.

A realistic workflow uses all three in sequence: write and test a pipeline on a laptop with a tiny subsampled FASTQ file (say, 100,000 reads), move it to HPC or cloud for the full run across all samples, and bring summary tables and plots back to the laptop for final figures. Module 15 (Workflow Managers and Reproducibility) covers how tools like Snakemake and Nextflow let you write a pipeline once and run it unchanged on all three.

1.2 The command line you actually need

Nearly every genomics file format is plain text, or a compressed/binary variant of plain text (Module 2 covers formats in detail). That makes the Unix shell — a text-based interface to the operating system — the natural tool for genomics: you are almost always filtering, counting, or reshaping lines of text, and the shell was designed for exactly that, decades before genomics existed.

Filesystem basics. pwd prints your current directory. ls -lh lists files with human-readable sizes. cd path changes directory (cd .. goes up one level, cd ~ goes home). mkdir -p project/data creates nested directories in one step. cp, mv, rm copy, move, and delete — rm has no undo and no trash can, so rm -rf some_dir deletes a directory tree permanently and immediately. Always ls a wildcard pattern before you rm it.

Pipes and why they fit genomics. A pipe (|) sends the output of one command directly into the input of the next, without ever writing an intermediate file to disk. This matters because genomics files are often tens of gigabytes; writing and re-reading intermediates repeatedly wastes disk space and time. A single pipeline can stream a 30 GB BAM file (a compressed alignment file, covered in Module 2) through several filtering and counting steps and never touch disk except for the final answer.

# Count how many reads in a FASTQ file are duplicates of the first 20 bases
# (a crude proxy for library complexity — not a real deduplication method)
zcat sample_R1.fastq.gz | awk 'NR%4==2' | cut -c1-20 | sort | uniq -c | sort -rn | head

Reading that pipeline left to right: zcat decompresses the gzipped FASTQ to standard output; awk 'NR%4==2' keeps only line 2 of every 4 (FASTQ records are 4 lines each, and line 2 is the sequence); cut -c1-20 keeps the first 20 characters of each sequence; sort groups identical strings together (required before uniq can count them); uniq -c counts consecutive identical lines; sort -rn sorts those counts numerically, largest first; head shows the top 10. Nothing here touches disk until you redirect it with >.

The core text tools, with what each one is actually for:

A worked combination: count variants per chromosome in a VCF file (variant call format, covered in Module 7), excluding header lines:

grep -v '^#' variants.vcf | cut -f1 | sort | uniq -c | sort -rn
#   152340 chr1
#   148221 chr2
#    ...

Loops, for when you need to repeat a command over many files:

for f in *.fastq.gz; do
    fastqc "$f" -o qc_reports/
done

Quoting "$f" matters — unquoted, a filename containing a space breaks the loop into two arguments.

Remote work: ssh, scp, rsync. ssh user@cluster.university.edu opens a remote shell on a cluster or cloud instance. scp localfile.fastq.gz user@cluster:/scratch/me/ copies a single file across, but re-copies the whole file even if 99% is unchanged. rsync -avz --progress data/ user@cluster:/scratch/me/data/ is almost always the better choice: it compares files and only transfers differences, resumes interrupted transfers, and preserves permissions and timestamps (-a for archive mode, -v for verbose, -z to compress in transit).

Keeping jobs alive after you log off. An ssh session dies the moment your connection drops, and by default any command you started dies with it. Two standard fixes: nohup long_command.sh > log.txt 2>&1 & detaches the command from the terminal and redirects both standard output and standard error (2>&1) to a log file, so it survives logout. tmux (or the older screen) is more flexible: it creates a persistent terminal session on the remote machine that you can detach from (Ctrl-b d in tmux) and reattach to later (tmux attach), even from a different physical computer, picking up exactly where you left off — including scrollback and any running interactive program.

Disk housekeeping. df -h shows free space per filesystem; du -sh */ shows how much space each subdirectory in the current folder uses, sorted by nothing by default — pipe through sort -rh to see the biggest first. Running out of disk mid-pipeline is one of the most common ways a long genomics job silently corrupts its own output, because many tools do not check write-return codes carefully.

Permissions. ls -l shows a string like -rw-r--r--: the first character is file type (- plain file, d directory), then three groups of rwx (read/write/execute) for owner, group, and everyone else. chmod 750 script.sh sets owner to read/write/execute, group to read/execute, others to nothing (the digits are a bitmask: read=4, write=2, execute=1, summed per group). On shared HPC systems, group permissions matter because lab members typically share a group, and institutional policy will restrict "others" access to protect patient or embargoed data.

1.3 Conda, mamba, and four complete environments

Why environment managers exist. Bioinformatics tools are written in different languages (C, C++, Python, Java, R), each with its own dependency versions, and different tools in the same pipeline often need incompatible versions of the same shared library. Conda solves this by creating isolated environments — self-contained sets of packages with their own binaries — so that gatk4 in one environment can depend on a different Java version than nextflow in another, with zero conflict. Mamba is a drop-in reimplementation of the conda command with a much faster dependency solver (minutes instead of sometimes hours for large environments); the file format and commands are otherwise identical, so everything below works with either conda env create or mamba env create.

Channel order and why it matters. A "channel" is a package repository. conda-forge is a large, general-purpose, community-maintained channel with up-to-date builds of core scientific Python/R packages. bioconda is a channel specifically for bioinformatics tools (aligners, variant callers, format utilities) that depends on packages from conda-forge for its own builds. Channel order in an environment.yml sets priority for resolving which channel to pull a package from when more than one channel provides it. The convention — conda-forge listed before bioconda, with defaults excluded entirely — exists because mixing in the defaults channel (Anaconda Inc.'s own channel) alongside bioconda is the single most common cause of unsolvable or silently broken environments in genomics; bioconda's own documentation mandates conda-forge > bioconda > nothing else.

The pip-inside-conda rule. Conda and pip (Python's native package installer) are two different package managers with no awareness of each other by default. If you pip install into a conda environment, pip can silently overwrite a library that conda installed and tracked, leaving conda's internal record of "what is installed" wrong — the next conda update or conda install may then break the environment because conda thinks a package is at one version when pip actually replaced it. The rule: install everything you can via conda/mamba first, and only use pip for packages genuinely unavailable on conda-forge/bioconda, listed under a pip: block at the end of the environment.yml so conda resolves and installs its own packages first, then hands the remainder to pip inside the already-built environment.

# environment-genomics.yml
name: genomics
channels:
  - conda-forge
  - bioconda
dependencies:
  - python=3.11
  - fastqc=0.12.1
  - fastp=0.23.4
  - bwa=0.7.17
  - minimap2=2.27
  - samtools=1.20
  - bcftools=1.20
  - gatk4=4.5.0.0
  - star=2.7.11b
  - salmon=1.10.2
  - subread=2.0.6
  - multiqc=1.21
  - snakemake=8.14.0
  - nextflow=24.04.2

# environment-scomics.yml
name: scomics
channels:
  - conda-forge
  - bioconda
dependencies:
  - python=3.10
  - numpy
  - scipy
  - pandas
  - scanpy=1.10.1
  - anndata=0.10.7
  - leidenalg=0.10.2
  - harmonypy=0.0.10
  - squidpy=1.4.1
  - igraph
  - pip
  - pip:
      - scvi-tools==1.1.2
      - scanorama==1.7.4
      - cell2location==0.1.4

# environment-pathology.yml
name: pathology
channels:
  - conda-forge
  - bioconda
dependencies:
  - python=3.10
  - openslide=4.0.0
  - openslide-python=1.3.1
  - tiffslide=2.4.0
  - opencv=4.9.0
  - scikit-image=0.23.1
  - numpy
  - pandas
  - matplotlib
  - pip
  - pip:
      - stardist==0.9.1
      - cellpose==3.0.10
      - timm==1.0.3
      - tensorflow==2.16.1

# environment-dlbio.yml
name: dlbio
channels:
  - conda-forge
  - bioconda
dependencies:
  - python=3.11
  - pytorch
  - pytorch-cuda=12.1
  - rdkit=2024.03.1
  - numpy
  - pandas
  - scikit-learn
  - pip
  - pip:
      - lightning==2.2.4
      - transformers==4.41.2
      - peft==0.11.1
      - deepchem==2.8.0

A few points apply to all four files. First, pinning exact versions (samtools=1.20 rather than a bare samtools) is deliberate: an unpinned install resolved a year from now can silently pick a newer version with changed default behaviour or a renamed flag, breaking a pipeline that worked when you wrote it — pin versions for anything you intend to publish or hand to a collaborator. Second, pytorch-cuda=12.1 in the dlbio file requests the CUDA-enabled build; on a machine with no NVIDIA GPU, drop that line and conda will install a CPU-only build instead, which still runs everything, just slower. Third, build each environment with mamba env create -f environment-genomics.yml (substitute the filename); update it later with mamba env update -f environment-genomics.yml --prune to remove packages dropped from the file as well as add new ones. Finally, keep one environment.yml per project checked into the project's version control (Module 16 covers Git), not one giant environment with everything installed — a narrow environment solves faster, conflicts less, and makes it obvious to a collaborator exactly what the analysis depends on.

1.4 Python for biologists

Python is a general-purpose programming language. In bioinformatics it is the default choice when you need to glue together file formats, write custom parsers, build machine learning pipelines, or scale an analysis beyond what a point-and-click tool offers. The value is not the language syntax itself — it is a small set of libraries that give you arrays, tables, and plots, plus domain libraries that understand biological file formats.

1.4.1 The four-library core

These four form a stack: pandas holds your data, numpy does the math under the hood, and matplotlib/seaborn draw it. You will import this combination in almost every analysis script you write.

1.4.2 Vectorisation: the single most important habit

A for loop over rows of a table, computed one scalar at a time, is the single most common performance mistake made by biologists coming from spreadsheet or script backgrounds. Vectorisation (applying one operation to an entire array at once, letting the library's compiled C code do the looping) is routinely 10-100x faster and is also shorter to write.

import numpy as np
import pandas as pd

counts = pd.DataFrame(
    np.random.poisson(lam=20, size=(5, 4)),
    index=[f"gene{i}" for i in range(5)],
    columns=["sample1", "sample2", "sample3", "sample4"],
)

# BAD: explicit Python loop, row by row
log_counts_slow = []
for i in range(counts.shape[0]):
    row = []
    for j in range(counts.shape[1]):
        row.append(np.log2(counts.iloc[i, j] + 1))
    log_counts_slow.append(row)

# GOOD: vectorised — one call, operates on the whole array
log_counts = np.log2(counts + 1)
# counts.shape == (5, 4); log_counts has identical shape and index/columns preserved

The vectorised version is not just faster — it is also easier to read, because it states the transformation once instead of describing how to iterate.

1.4.3 groupby-apply: the workhorse of tabular analysis

Most biological questions reduce to "compute something, per group." Per-gene mean expression across replicates. Per-sample total read count. Per-patient maximum variant allele frequency. pandas' groupby implements the split-apply-combine pattern: split the table into groups by a key, apply a function to each group, combine the results back into one table.

import pandas as pd

df = pd.DataFrame({
    "gene": ["BRCA1", "BRCA1", "TP53", "TP53", "EGFR", "EGFR"],
    "sample": ["s1", "s2", "s1", "s2", "s1", "s2"],
    "tpm": [12.4, 15.1, 102.3, 98.7, 5.0, 4.8],
})

# mean TPM per gene, across samples
gene_means = df.groupby("gene")["tpm"].mean()
# gene_means:
# gene
# BRCA1    13.75
# EGFR      4.90
# TP53    100.50
# Name: tpm, dtype: float64

# custom per-group function via apply
def coefficient_of_variation(x):
    return x.std() / x.mean()

cv_per_gene = df.groupby("gene")["tpm"].apply(coefficient_of_variation)

.apply is flexible but slower than built-in aggregations (.mean(), .sum(), .agg()) because it calls a Python function per group instead of running a vectorised routine. Use .agg() with built-in names or multiple functions at once when you can:

summary = df.groupby("gene")["tpm"].agg(["mean", "std", "min", "max"])

1.4.4 Tidy vs wide: the shape decision that determines everything downstream

This is the most consequential formatting decision in a tabular analysis, and it recurs in every module that touches counts matrices, metadata, or statistical tests.

• Wide format: one row per feature (gene), one column per sample. This is how expression matrices are naturally produced by quantification tools (featureCounts, salmon, Cell Ranger). Good for matrix algebra (PCA, clustering, linear models expressed as $X\beta$).
• Tidy (long) format: one row per observation — a (gene, sample, value) triple, with additional columns for metadata (condition, batch, timepoint). Good for groupby, for joining with sample metadata, and for seaborn/ggplot2 plotting, both of which expect one row per data point.

Converting between the two is routine and you should be fluent in both directions:

# wide -> long ("melt")
wide = pd.DataFrame({
    "gene": ["BRCA1", "TP53"],
    "s1": [12.4, 102.3],
    "s2": [15.1, 98.7],
})
long = wide.melt(id_vars="gene", var_name="sample", value_name="tpm")
#    gene sample    tpm
# 0 BRCA1     s1   12.4
# 1  TP53     s1  102.3
# 2 BRCA1     s2   15.1
# 3  TP53     s2   98.7

# long -> wide ("pivot")
back_to_wide = long.pivot(index="gene", columns="sample", values="tpm")

The rule of thumb used throughout this course: store and compute in wide when the operation is matrix algebra (normalisation, PCA, differential expression model fitting); convert to tidy/long for plotting, filtering by metadata, and ad hoc summaries. Module 6 (RNA-seq Quantification and Differential Expression) and Module 11 (Statistical Foundations) depend on you being able to move fluently between these two shapes without re-deriving the logic each time.

1.4.5 Biopython quickstart

Biopython is the standard library for parsing biological sequence file formats, talking to NCBI, and doing simple sequence manipulations in Python. You will reach for it to read FASTA/GenBank files, translate sequences, or compute reverse complements — not for high-performance alignment, which belongs to dedicated compiled tools (Module 2, Alignment).

from Bio import SeqIO
from Bio.Seq import Seq

# iterate over a FASTA file without loading it all into memory
for record in SeqIO.parse("genome.fasta", "fasta"):
    print(record.id, len(record.seq))
    # NC_000913.3 4641652

# basic sequence operations
s = Seq("ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG")
print(s.reverse_complement())
# CTATCGGGCACCCTTTCAGCGGCCCATTACAATGGCCAT
print(s.translate(to_stop=True))
# MAIVMGR*  (translation runs to stop codon; '*' only shown if to_stop=False)

# write filtered records back out
long_records = (r for r in SeqIO.parse("genome.fasta", "fasta") if len(r.seq) > 1000)
SeqIO.write(long_records, "genome_long_only.fasta", "fasta")

SeqIO.parse returns an iterator, not a list — for large genome or multi-FASTA files (millions of records, e.g. a transcriptome), iterate rather than calling list(...) on it, or you will exhaust memory.

1.4.6 pysam quickstart

pysam is a Python wrapper around the htslib C library (the same engine behind samtools, bcftools, tabix). It lets you read and query BAM/SAM (aligned reads), VCF/BCF (variants), and indexed FASTA files directly from Python, without shelling out to samtools for every operation.

import pysam

# reading an indexed BAM file (requires a .bai index alongside it)
bam = pysam.AlignmentFile("sample.bam", "rb")

# fetch reads overlapping a region (0-based, half-open, matches BAM convention)
for read in bam.fetch("chr1", 1000000, 1000100):
    print(read.query_name, read.reference_start, read.mapping_quality)

# compute per-base coverage over a region
coverage = bam.count_coverage("chr1", 1000000, 1000100)
# coverage is a tuple of 4 arrays (A, C, G, T counts), each length 100

bam.close()

# reading a VCF
vcf = pysam.VariantFile("variants.vcf.gz")
for rec in vcf.fetch("chr1", 1000000, 2000000):
    print(rec.chrom, rec.pos, rec.ref, rec.alts, rec.qual)

Two details that cause silent bugs: BAM/BED coordinates are 0-based half-open (fetch(start, end) includes start, excludes end), while VCF and most human-readable genome browsers use 1-based inclusive coordinates. pysam exposes rec.pos in VCF's native 1-based convention but read.reference_start in BAM's 0-based convention — mixing these up when writing coordinates between formats is one of the most common sources of off-by-one errors in genomics code (Module 2 covers coordinate systems in depth).

1.5 R for biologists

R is a statistical programming language developed specifically for data analysis, and it has a second, independent ecosystem beyond CRAN (the general package repository) called Bioconductor, built specifically for biological data: expression matrices with gene/sample annotations attached, genomic ranges, variant data, and standardised statistical methods for them (differential expression, batch correction, pathway enrichment). This is why R remains dominant in parts of genomics even though Python has caught up or overtaken it elsewhere.

1.5.1 tidyverse: the data manipulation grammar

The tidyverse is a coordinated set of R packages (dplyr, tidyr, ggplot2, readr, purrr, among others) built around one shared idea: tidy data, where every row is an observation and every column is a variable — the same tidy/long concept introduced in 1.4.4, developed first and most thoroughly in the R ecosystem.

library(tidyverse)

df <- tibble(
  gene   = c("BRCA1", "BRCA1", "TP53", "TP53"),
  sample = c("s1", "s2", "s1", "s2"),
  tpm    = c(12.4, 15.1, 102.3, 98.7)
)

# the pipe (%>% or, in base R >= 4.1, |>) chains operations left to right
summary_tbl <- df %>%
  group_by(gene) %>%
  summarise(mean_tpm = mean(tpm), sd_tpm = sd(tpm))
# # A tibble: 2 x 3
#   gene  mean_tpm sd_tpm
#   <chr>    <dbl>  <dbl>
# 1 BRCA1    13.75   1.91
# 2 TP53    100.50   2.55

wide <- df %>% pivot_wider(names_from = sample, values_from = tpm)
long_again <- wide %>% pivot_longer(cols = -gene, names_to = "sample", values_to = "tpm")

pivot_wider/pivot_longer are the direct R analogues of pivot/melt in pandas — the same wide-vs-tidy decision from 1.4.4 applies verbatim in R.

1.5.2 Bioconductor and BiocManager

Bioconductor packages are not installed from CRAN. They are installed through BiocManager, a CRAN package whose job is to install Bioconductor packages at versions that are mutually compatible (Bioconductor releases packages in coordinated twice-yearly bundles, unlike CRAN's rolling releases).

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install(c("SummarizedExperiment", "DESeq2", "limma"))
BiocManager::valid()  # checks that installed Bioconductor packages match your R/Bioconductor version

1.5.3 SummarizedExperiment and SingleCellExperiment: data containers that carry their own metadata

A plain matrix of counts (genes x samples) loses information the moment you separate it from the sample metadata (condition, batch, sequencing depth) and feature metadata (gene length, chromosome). Bioconductor's core data structure, SummarizedExperiment, bundles all three together in one object, so that subsetting the matrix automatically subsets the matching metadata rows — this is the Bioconductor analogue of the pandas DataFrame with an index, but extended to carry two metadata tables at once (one for rows/features, one for columns/samples).

library(SummarizedExperiment)

counts_mat <- matrix(rpois(20, lambda = 20), nrow = 5, ncol = 4,
                      dimnames = list(paste0("gene", 1:5), paste0("sample", 1:4)))

col_data <- DataFrame(condition = c("control", "control", "treated", "treated"),
                       batch     = c("A", "B", "A", "B"))
rownames(col_data) <- colnames(counts_mat)

se <- SummarizedExperiment(assays = list(counts = counts_mat), colData = col_data)

se                      # shows dimensions, assay names, colData columns
assay(se, "counts")     # the matrix itself
colData(se)             # sample metadata, row-aligned to columns of the matrix
se[, se$condition == "treated"]   # subsetting keeps matrix and metadata in sync

SingleCellExperiment (from the SingleCellExperiment package) extends SummarizedExperiment for single-cell data: it adds slots for reduced dimensions (PCA, UMAP coordinates) and alternative experiments (e.g. spike-in or antibody-derived tag counts) alongside the main counts matrix. Module 8 (Single-Cell Genomics) uses SingleCellExperiment objects (or Python's AnnData, the equivalent container in the scanpy ecosystem) as the standard unit of data throughout — knowing this container now means Module 8 can focus on the biology and statistics instead of re-explaining the data structure.

1.5.4 ggplot2: the grammar of graphics

ggplot2 builds plots by layering independent components onto a shared coordinate system, following a model called the "grammar of graphics": you declare data, an aesthetic mapping (which column maps to x, y, colour, shape), and one or more geometric layers (points, lines, boxplots), and the library composes them.

library(ggplot2)

ggplot(df, aes(x = gene, y = tpm, fill = sample)) +
  geom_boxplot() +
  labs(title = "Expression by gene", y = "TPM") +
  theme_minimal()

The seaborn equivalent in 1.4 covers the same statistical plot types; the practical difference is that ggplot2's layered + syntax makes it easy to add or swap one component (a different geom, a facet, a theme) without rewriting the whole call, which is why R retains an edge for publication-quality exploratory plotting even among groups that do their main analysis in Python.

1.5.5 When to choose R over Python in this field

The honest answer for most labs is: learn both at a working level. Bulk RNA-seq differential expression and most classical Bioconductor statistical methods are easiest in R; everything involving general scripting, large-scale data engineering, or machine learning is easiest in Python. Module 6 and Module 11 will give you both an R and a Python path for the core statistical workflows so you are not locked into one ecosystem.

1.6 Jupyter, VS Code, RStudio, Quarto, and notebook hygiene

These are the four environments you will actually type code into. They solve different problems and are not mutually exclusive — many bioinformaticians use RStudio for R analysis, VS Code for pipeline/script development, and Quarto to produce the final shared report, within the same project.

1.6.1 Notebook hygiene

Notebooks (Jupyter or RStudio's R Notebooks/Quarto documents in interactive mode) are excellent for exploration and terrible for reproducibility if used carelessly, because the displayed output reflects the order cells were run, not the order they appear on the page. A cell you edited and reran out of sequence can leave the notebook showing results that no longer match the code above them.

Rules that prevent this, in order of importance:

1. Before trusting or sharing a notebook, restart the kernel and run all cells top to bottom (Jupyter: "Restart & Run All"; RStudio: "Restart R" then source the document). If it does not reproduce the same output cleanly, the notebook is lying about its own state.
2. Number or name cells by logical step, not by edit history — do not leave three abandoned attempts above the one that worked.
3. Do not commit notebooks with large embedded outputs (plots, big printed tables) to git unversioned; .ipynb files store output as embedded base64/JSON, which bloats repositories and produces unreadable diffs. Use nbstripout (strips output before commit) or move final, stable analyses into plain .py/.R scripts or Quarto documents.
4. One notebook, one question. A notebook that starts as "QC the counts matrix" and drifts into "also try five different normalisation methods and a clustering experiment" becomes unreadable and unrerunnable. Split it.
5. Pin your environment (Module 1.2/1.3's conda/renv discipline applies to notebooks too) — a notebook that ran six months ago under a different package version is not guaranteed to run the same way today.
6. For anything meant to be read by someone else — a collaborator, a reviewer, your future self in a year — render to Quarto or a script, not a live notebook. A .qmd file renders deterministically from a clean environment every time, which is exactly the guarantee an interactive notebook does not give you by default.

1.7 Git and GitHub essentials for a scientist

Git is a version control system (a program that records snapshots of a folder's contents over time, so you can go back, compare, and merge changes). GitHub is a website that hosts Git repositories (a repository, or "repo," is a folder tracked by Git) and adds collaboration features: issues, pull requests, code review, and a web UI. You do not need to be a software engineer to use either. You need about fifteen commands and one mental model.

The mental model. Git keeps three areas: the working directory (your files as they sit on disk), the staging area (a holding pen for changes you are about to record), and the history (a chain of snapshots called commits). You edit files, git add the ones you want recorded, then git commit to seal a snapshot with a message. git push sends your local history to a remote copy (GitHub); git pull fetches and merges someone else's new commits into yours.

# one-time setup on a new machine
git config --global user.name  "Jane Researcher"
git config --global user.email "jane@example.edu"
git config --global init.defaultBranch main

# start tracking a new project
mkdir rnaseq-liver-project && cd rnaseq-liver-project
git init
echo "results/\n*.bam\n*.fastq.gz\n.Rhistory\n__pycache__/" > .gitignore
git add .gitignore README.md
git commit -m "Initial project skeleton"

# connect to GitHub (create the empty repo on github.com first)
git remote add origin git@github.com:janeresearcher/rnaseq-liver-project.git
git push -u origin main

A .gitignore file tells Git which files never to track — raw data, large intermediate files (BAM, FASTQ), credentials, and editor junk. Never commit raw sequencing data or credentials to Git. Git repositories are meant for code, small text files, and configuration, not multi-gigabyte binaries; GitHub rejects files over 100 MB outright and degrades badly well before that. For large files that must be versioned, use Git LFS (Large File Storage, an extension that stores a pointer in Git and the actual bytes elsewhere) or, better, keep data out of Git entirely and track it with a data-manifest (a checksummed list of file paths and hashes) instead.

Daily workflow.

git status                       # what changed since the last commit?
git add scripts/qc_filter.py     # stage one file
git add -p                       # stage changes interactively, chunk by chunk
git commit -m "Add adapter-trimming step to QC script"
git log --oneline -10            # last 10 commits, one line each
git diff HEAD~1                  # what changed in the last commit?

Branches let you work on a change in isolation without disturbing the working version. A branch is a movable pointer to a commit; creating one is instant and cheap.

git checkout -b feature/add-deseq2-contrast
# ... edit files, commit ...
git push -u origin feature/add-deseq2-contrast
# open a Pull Request (PR) on github.com: a request to merge this branch into main,
# with a diff view and a place for a collaborator to comment before it's merged

A pull request is GitHub's review unit: it shows the diff, lets co-authors comment line-by-line, and runs automated checks (continuous integration, CI) before the branch is merged into main. For a solo analyst, PRs are still useful as a changelog and a place to leave yourself notes on why a change was made, not just what changed.

Handling a merge conflict (it will happen): Git marks the clashing lines with <<<<<<<, =======, >>>>>>> inside the file. Open the file, decide which version (or combination) is correct, delete the markers, then git add the file and git commit to close the merge. There is no magic command that picks correctly for you — that judgment is the point of version control.

Reproducibility payoff. Every analysis script, notebook, and pipeline definition should live in Git. When you report a result, you should be able to say "this figure was produced by commit a3f9c12 of this repository" — a precise, checkable claim that "I ran some R code last Tuesday" is not. Tag the exact commit used for a paper's analyses:

git tag -a v1.0-submission -m "State of the analysis at manuscript submission"
git push origin v1.0-submission

Recommended project directory layout

A consistent layout means you (and collaborators, and Reviewer 2 eighteen months later) can find anything in ten seconds without a README scavenger hunt. This is a convention, not a law — adapt it, but adopt some convention and stick to it across projects.

rnaseq-liver-project/
├── README.md                # what this project is, how to reproduce it
├── .gitignore
├── environment.yml          # conda/mamba spec, or renv.lock / requirements.txt
├── config/
│   └── samples.tsv          # sample sheet: sample_id, condition, fastq paths
├── data/
│   ├── raw/                 # untouched inputs — read-only, never edited (often not in Git)
│   └── external/            # reference genomes, annotation, downloaded once
├── scripts/                 # or src/ — one script per pipeline step, numbered
│   ├── 01_qc_trim.sh
│   ├── 02_align.sh
│   ├── 03_count.sh
│   └── 04_deseq2.R
├── notebooks/                # exploratory Jupyter/RMarkdown — not the final pipeline
├── results/                  # generated outputs: tables, BAMs, counts (not in Git)
├── figures/                  # generated plots (not in Git, or tracked with LFS)
└── docs/                      # manuscript drafts, method notes

Rules of thumb that save real pain later: raw data is read-only and never edited in place; every generated file can be deleted and regenerated from scripts/ plus data/; file and sample identifiers in config/samples.tsv are the single source of truth that every script reads from, so renaming a sample means editing one file, not twenty.

1.8 HPC: SLURM basics, resource estimation, and job arrays

A high-performance computing (HPC) cluster is a set of networked computers (nodes) sharing storage, managed by a scheduler that decides whose job runs where and when. You do not get to run code directly on a cluster's compute nodes by typing commands interactively (usually) — you log into a login node (for editing files and submitting jobs, not for heavy computation), and you submit jobs (descriptions of a program to run plus the resources it needs) to a queue. SLURM (Simple Linux Utility for Resource Management) is the most common scheduler in academic HPC; others include PBS/Torque and LSF, with broadly similar concepts under different command names.

Core SLURM commands.

A worked sbatch script for aligning one RNA-seq sample with STAR:

#!/bin/bash
#SBATCH --job-name=star_align
#SBATCH --partition=general
#SBATCH --nodes=1
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=8
#SBATCH --mem=32G
#SBATCH --time=02:00:00
#SBATCH --output=logs/star_%j.out        # %j = job ID
#SBATCH --error=logs/star_%j.err
#SBATCH --mail-type=END,FAIL
#SBATCH --mail-user=jane@example.edu

set -euo pipefail
module load star/2.7.11a               # cluster-specific module system

STAR \
  --runThreadN "${SLURM_CPUS_PER_TASK}" \
  --genomeDir ref/star_index_GRCh38 \
  --readFilesIn data/raw/sampleA_R1.fastq.gz data/raw/sampleA_R2.fastq.gz \
  --readFilesCommand zcat \
  --outSAMtype BAM SortedByCoordinate \
  --outFileNamePrefix results/sampleA_

Submit with sbatch align_sampleA.sh. The directives after #SBATCH are read by the scheduler, not the shell; they must appear before any executable line. --cpus-per-task=8 requests eight CPU threads on one node for one task; --mem=32G requests total memory for the job (not per-CPU, unless you use --mem-per-cpu); --time=02:00:00 is a hard wall-clock limit — exceed it and SLURM kills the job, no matter how close to finishing it was.

Resource estimation. Guessing resources is the single biggest source of wasted queue time: ask for too little memory and the job is killed (OutOfMemory, oom-kill); ask for too much and the scheduler makes you wait far longer for a node that satisfies the request, and you burn allocation (many HPC centers bill or quota by cpu-hours × requested-memory, not actual usage). The right workflow is empirical: run one representative sample with generous resources, check actual usage with sacct, then size the rest of the batch accordingly.

sacct -j 184231 --format=JobID,Elapsed,MaxRSS,ReqMem,State
# JobID       Elapsed   MaxRSS     ReqMem   State
# 184231      00:41:12  28.3G      32G      COMPLETED

MaxRSS (maximum resident set size — peak real memory used) of 28.3 GB against a 32 GB request is a sane margin (roughly 15-20% headroom). Rough starting points for common genomics steps, to be corrected against your own sacct numbers:

Job arrays run the same script over many inputs (samples) as one logical submission, each with a different index available as $SLURM_ARRAY_TASK_ID. This is the standard pattern for "do this to all 96 samples" instead of writing 96 scripts.

#!/bin/bash
#SBATCH --job-name=fastp_array
#SBATCH --array=1-96%12          # 96 tasks, at most 12 running at once
#SBATCH --cpus-per-task=4
#SBATCH --mem=4G
#SBATCH --time=00:30:00
#SBATCH --output=logs/fastp_%A_%a.out   # %A = array job ID, %a = task index

sample=$(sed -n "${SLURM_ARRAY_TASK_ID}p" config/samples.tsv | cut -f1)
fastp -i data/raw/${sample}_R1.fastq.gz -I data/raw/${sample}_R2.fastq.gz \
      -o results/${sample}_R1.trim.fastq.gz -O results/${sample}_R2.trim.fastq.gz \
      --thread 4

%12 caps concurrency so you don't monopolize the partition or hammer shared filesystem I/O. Each array element is accounted and scheduled independently, so sacct -j <arrayid> shows every task's resource usage separately — invaluable for spotting the one sample that used 10x the memory of its siblings (often a sign of a corrupted file or an outlier library size, not a resource problem to just throw more RAM at).

The free-tier reality

Not every reader has an HPC allocation. A realistic hierarchy of free and low-cost compute, with honest limits:

A pragmatic strategy: prototype on Colab or a laptop with a down-sampled dataset, move to institutional HPC or Galaxy for the full production run, and reserve commercial cloud credits for the moments you need elastic scale (a one-off 500-sample reanalysis, a GPU-hungry deep learning sweep covered in Module 9) rather than as a default.

1.9 Where public data lives, and how to get it

Public repositories are the backbone of reproducible genomics: raw reads go to sequence archives, processed expression data to expression archives, and clinical/cohort data to controlled-access archives with a review process. Knowing which repository holds which kind of thing, and which tool pulls it down, saves hours of searching.

sra-tools, the standard way to pull raw reads from SRA:

# download the SRA-format archive for one run accession
prefetch SRR12345678

# convert to gzipped paired FASTQ, splitting read 1 / read 2
fasterq-dump SRR12345678 --split-files --threads 8 -O fastq/
gzip fastq/SRR12345678_1.fastq fastq/SRR12345678_2.fastq

ffq (a lightweight metadata-and-URL fetcher) resolves an accession — SRA run, GEO series, ENA study — straight to downloadable URLs and JSON metadata, which is often faster than SRA's own tools because it points you at ENA's direct FASTQ mirrors (no conversion step needed):

pip install ffq
ffq --ftp SRR12345678 > urls.json
# parse urls.json for "url" fields, then:
wget -c $(python -c "import json;print(json.load(open('urls.json'))[0]['url'])")

Bulk ENA downloads by direct FTP pattern (ENA lays out paths predictably by accession):

# ENA FASTQ path pattern: era-fasta/vol1/fastq/<first6>/<00+lastdigit if 10+ digits>/<accession>/
wget -c "ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR123/078/SRR12345678/SRR12345678_1.fastq.gz"

GDC controlled data needs an authentication token downloaded from the GDC portal after your data access request is approved:

gdc-client download -m manifest.tsv -t gdc-user-token.txt -d gdc_downloads/

Decision table: "I want to do X → use this stack"

1.10 Common pitfalls and how to avoid them

1.11 Exercises

1. (Warm-up) Initialize a Git repository for a toy project, create a .gitignore that excludes *.bam and results/, make two commits, and push to a new GitHub repo. Deliverable: the GitHub URL and output of git log --oneline.
2. (Warm-up) Write an sbatch script that requests 2 CPUs, 4 GB memory, and 10 minutes, and runs echo "hello from $(hostname)". Submit it and retrieve the output file. Deliverable: the script and the contents of the .out file.
3. (Core) Take a GEO series accession (pick any public one from a paper you know) and trace it to its linked SRA/BioProject accession, then write the ffq or ENA FTP command you would use to fetch raw FASTQ for one run. Deliverable: the accession chain (GSEnnnn → GSMnnnn → SRRnnnn) and the exact download command.
4. (Core) Convert exercise 2's single job into a job array of 5 tasks, each printing its own $SLURM_ARRAY_TASK_ID to a separate log file, capped at 2 concurrent tasks. Deliverable: the modified script and a listing of the 5 output files.
5. (Core) For a 16-sample human RNA-seq experiment, estimate total CPU-hours and total GB-hours needed for STAR alignment, using the per-sample figures in the resource table in this module. State your assumptions. Deliverable: a short table and the arithmetic.
6. (Stretch) Design a Git branching workflow for a two-person analysis project: one person writes preprocessing scripts, the other writes statistical analysis. Specify branch names, when to merge, and how you'd structure pull requests to catch a bug before it reaches main. Deliverable: a one-page written plan.
7. (Stretch) You need individual-level genotype data from a European cohort study for a project with institutional ethics approval. Outline, step by step, the real-world process from "I found the study" to "I have the data on a secure machine," naming the specific repository and access mechanism. Deliverable: a numbered procedure referencing EGA.

Solutions / hints

1. git init, then printf '*.bam\nresults/\n' > .gitignore, git add .gitignore && git commit -m "Add gitignore", add a file and git commit -m "Add script", then git remote add origin <url> && git push -u origin main. git log --oneline should show two commits.
2. Script needs #SBATCH --cpus-per-task=2, #SBATCH --mem=4G, #SBATCH --time=00:10:00, and the body echo "hello from $(hostname)". The .out file (named via #SBATCH --output=hello_%j.out) contains one line: hello from <nodename>.
3. Example chain pattern: a GEO Series GSE123456 page lists Samples GSM3456789 etc.; each GSM page links "SRA" to a Run accession like SRRxxxxxxx under a BioProject PRJNAxxxxxx. Command: ffq --ftp SRRxxxxxxx > urls.json then wget the FASTQ URL inside it, or the direct ENA FTP path using the SRR's first-6-characters/padding convention shown in section 1.9.
4. Replace single-job directives with #SBATCH --array=1-5%2, change #SBATCH --output=hello_%A_%a.out, and the body becomes echo "hello from task ${SLURM_ARRAY_TASK_ID} on $(hostname)". Five files hello_<jobid>_1.out through _5.out should appear, at most 2 running simultaneously (visible via squeue while running).
5. Using 8-12 CPUs and 32 GB per sample at ~30 min wall time: CPU-hours ≈ 16 samples × 10 CPUs × 0.5 h = 80 CPU-hours; GB-hours ≈ 16 × 32 GB × 0.5 h = 256 GB-hours. State explicitly that this assumes uniform sample size and that real libraries vary, so a 20-30% safety margin on time and memory is reasonable.
6. A reasonable plan: main always reflects a working, reproducible state; each person works on feature/preprocessing and feature/deseq2-analysis branches; preprocessing is merged first via a PR that the analysis-author reviews (checking output file formats match what the analysis script expects); the analysis branch is rebased on updated main before its own PR; nothing is merged into main without at least one read-through from the other person, catching format mismatches before they become silent bugs downstream.
7. (1) Identify the study and its EGA accession (e.g. EGASxxxxxxx) via its publication or the EGA website search. (2) Register an EGA account. (3) Submit a Data Access Request to the study's Data Access Committee (DAC), including your institutional ethics approval and a description of intended use. (4) Wait for DAC approval, which may take weeks and may be time-limited or project-scoped. (5) On approval, receive credentials for pyEGA3. (6) Download onto a secure, access-controlled machine that meets your institution's data security requirements — not a laptop or shared compute node — per the terms of the data access agreement.

1.12 Key takeaways

• Git records snapshots (commits) of code; GitHub hosts and adds review (pull requests) on top.
• Never commit raw data or credentials to Git; use .gitignore from the first commit.
• A consistent project directory layout (raw data read-only, scripts numbered, generated output regenerable) saves far more time than it costs to set up.
• SLURM jobs are described in #SBATCH directives in a script and submitted with sbatch; squeue shows status, sacct shows what resources a finished job actually used.
• Resource requests should be set empirically from sacct output on a representative run, not guessed.
• Job arrays (--array=1-N%K) are the standard way to run the same script over many samples with controlled concurrency.
• Colab, Galaxy, and institutional HPC are free or near-free; Terra and commercial cloud generally require funded credits and bill per use.
• SRA/ENA hold raw reads; GEO/ArrayExpress hold mostly processed expression data; know which one a given accession points to before choosing a download tool.
• sra-tools (prefetch/fasterq-dump) and ffq are the standard command-line routes to raw sequencing data; ENA's predictable FTP paths are often the fastest route.
• Controlled-access human data (UK Biobank, dbGaP, EGA) requires a formal access application and approval before any download is possible — there is no shortcut.
• Match the repository to the data type and the access tier before writing any download code.

1.13 Further reading

• Git: official documentation at git-scm.com; Pro Git book (Scott Chacon and Ben Straub), freely available online.
• GitHub: docs.github.com, particularly the guides on pull requests and .gitignore templates.
• SLURM: official documentation at slurm.schedmd.com, especially the sbatch, squeue, and job array pages.
• Google Colab: official FAQ and documentation at research.google.com/colaboratory.
• Galaxy Project: official documentation and training material at galaxyproject.org and training.galaxyproject.org.
• Terra: documentation at support.terra.bio.
• NCBI SRA and sra-tools: documentation and GitHub repository at github.com/ncbi/sra-tools.
• ENA: documentation at ebi.ac.uk/ena, including the FASTQ file path conventions.
• ffq: Pachter Lab GitHub repository, github.com/pachterlab/ffq.
• NCBI GEO: official documentation at ncbi.nlm.nih.gov/geo; Bioconductor GEOquery package documentation.
• EBI ArrayExpress / BioStudies: official documentation at ebi.ac.uk/biostudies.
• NCI Genomic Data Commons (GDC): docs.gdc.cancer.gov, including the gdc-client user guide.
• Human Cell Atlas Data Portal: data.humancellatlas.org.
• 10x Genomics public datasets: support.10xgenomics.com/datasets.
• ENCODE Project Portal: encodeproject.org, including its file download API documentation.
• UK Biobank Access Management System: ams.ukbiobank.ac.uk.
• dbGaP: official documentation at ncbi.nlm.nih.gov/gap, including data access request procedures.
• European Genome-phenome Archive (EGA): ega-archive.org, and the pyEGA3 client on GitHub.

Module 2 — Molecular Biology as a Data Model

In one paragraph. Every file you will ever load in a bioinformatics pipeline is a stylized description of a physical molecule — a strand of DNA, a spliced transcript, a folded protein — and the columns, coordinates, and flags in that file exist because of specific chemistry and cell biology. This module builds that biology from first principles: the chemistry of DNA, how genomes are packaged and expressed, how RNA is processed into the transcripts you sequence, and how proteins fold and get modified. The goal is that when you later see a BAM flag, a GTF exon, or a VCF genotype, you know exactly what real-world object it is standing in for.

Prerequisites: Module 1 (course orientation and computing setup). No prior biology assumed. You will be able to: Explain the chemical logic of base pairing and strand directionality; distinguish genome, chromosome, ploidy, and karyotype; state the central dogma and name its major exceptions; describe how a promoter, enhancer, and TAD differ and why distance in base pairs is not distance in regulatory influence; explain chromatin marks and DNA methylation as information layers on top of sequence; trace a pre-mRNA from capping through splicing to polyadenylation and explain why isoforms exist; justify why library prep protocols must choose which RNA species to capture; translate a coding sequence using the genetic code and explain codon usage bias; describe the four levels of protein structure and name common post-translational modifications. Time: 4-6 hours.

Figure 2.1 — From sample to claim. Every omics experiment is this pipeline. Each arrow discards information; the red labels are where it is lost or distorted. Nothing downstream can recover what an earlier stage threw away.

2.1 DNA chemistry and the double helix

DNA (deoxyribonucleic acid) is a polymer of four building blocks called nucleotides, each made of a deoxyribose sugar, a phosphate group, and one of four nitrogenous bases: adenine (A), thymine (T), guanine (G), cytosine (C). The sugars and phosphates alternate to form a backbone; the bases stick out sideways and pair with bases on a second, complementary strand through hydrogen bonds — A with T (two hydrogen bonds), G with C (three hydrogen bonds). Because G-C pairs have one more hydrogen bond, GC-rich DNA is more thermally stable, which matters directly for PCR primer design and for why GC-rich genomic regions are harder to sequence and assemble.

The two strands run in opposite chemical orientations (antiparallel). Each strand has a 5' end (named for the fifth carbon of the sugar, where a free phosphate sits) and a 3' end (the third carbon, where a free hydroxyl sits). Polymerases — the enzymes that build new DNA or RNA — can only add nucleotides to a 3' end, so synthesis always proceeds 5' to 3'. This single fact explains why one strand (the leading strand) replicates continuously and the other (the lagging strand) replicates in fragments, and it is why every sequence you will ever write in a file is implicitly read 5' to 3' unless stated otherwise.

The two strands wind around each other into a double helix, roughly 2 nm wide, making one full turn about every 10.5 base pairs (bp). This geometry is not decorative: the width of the major and minor grooves it creates is what allows proteins (transcription factors, polymerases) to read the base sequence without unwinding the helix, by inserting side chains into the grooves and sensing the pattern of hydrogen-bond donors and acceptors.

Failure mode: forgetting strandedness. A sequence and its reverse complement look unrelated as strings but describe the same physical DNA. Tools silently fail or give wrong answers when a BED interval's strand field is ignored, or when a primer is designed against the wrong strand.

2.2 Chromosomes, ploidy, karyotype

A chromosome is a single, very long DNA molecule (in humans, tens to hundreds of millions of bp) wound around proteins for packaging. The full set of DNA in a cell is its genome. Humans are diploid: somatic cells carry two copies of each autosome (chromosomes 1-22) plus two sex chromosomes (XX or XY), for 46 chromosomes total — one set inherited from each parent. Germ cells (sperm, egg) are haploid, carrying one copy of each chromosome, so that fertilization restores diploidy. Ploidy is not fixed across biology: many plants are polyploid (wheat is hexaploid), and some human tissues (liver, placenta) are naturally polyploid or aneuploid (abnormal chromosome number) as part of normal development — aneuploidy is only pathological in certain contexts, most famously trisomy 21 (Down syndrome, three copies of chromosome 21).

A karyotype is a visual, microscope-based inventory of an individual's chromosomes by number, size, and banding pattern, historically obtained by staining metaphase chromosomes (Giemsa staining, hence "G-banding"). It is the oldest cytogenetic data type and still the clinical standard for detecting large structural abnormalities (whole-chromosome gains/losses, large translocations) that short-read sequencing can miss or over-interpret. Karyotype notation (e.g., 46,XY,t(9;22)(q34;q11) for the BCR-ABL translocation in chronic myeloid leukemia) is a compact text format you will encounter in clinical genetics data even though it predates digital sequencing entirely.

2.3 The central dogma and its exceptions

The central dogma, as originally framed by Francis Crick, describes the directional flow of sequence information: DNA is transcribed into RNA, and RNA is translated into protein. As a default mental model this is correct and useful, but several well-established exceptions matter for interpreting real data:

Failure mode: treating "the genome" as the complete explanation for a phenotype. Sequence is necessary but not sufficient; which genes are switched on, in which cells, at what level, is a separate and equally large layer of information — the subject of the rest of this module.

2.4 Transcription and its machinery

Transcription is the synthesis of RNA from a DNA template, carried out by RNA polymerase. Humans have three main RNA polymerases with different jobs: Pol I makes ribosomal RNA (rRNA), Pol II makes messenger RNA (mRNA) and many regulatory RNAs, Pol III makes transfer RNA (tRNA) and other small RNAs. Pol II is the one most bioinformatics work centers on, because it produces the mRNA that becomes protein and most of what RNA-seq measures.

Pol II does not bind DNA alone. It requires general transcription factors (GTFs — TFIIA, TFIIB, TFIID, TFIIE, TFIIF, TFIIH) that assemble at the promoter into a pre-initiation complex. TFIID recognizes core promoter elements such as the TATA box (a TA-rich sequence about 25-30 bp upstream of the transcription start site in a subset of genes) and positions Pol II precisely at the transcription start site (TSS) — the genomic coordinate recorded as position 1 of a transcript's exon structure in annotation files. Sequence-specific transcription factors (proteins that bind particular DNA motifs, e.g., "GATA1 binds the motif GATA") then recruit or repress this machinery at specific genes, which is the molecular basis of every ChIP-seq and motif-enrichment analysis in later modules.

# Reverse-complementing a promoter motif to search both strands
# (transcription factors can bind either strand of the double helix)
def revcomp(seq: str) -> str:
    comp = str.maketrans("ACGTacgt", "TGCAtgca")
    return seq.translate(comp)[::-1]

motif = "TATAAA"          # canonical TATA box
print(revcomp(motif))      # -> 'TTTATA'

2.5 Promoters, enhancers, silencers, insulators, TADs, and 3D genome organization

A promoter sits immediately upstream of a gene's TSS and is where the general transcription machinery assembles. Enhancers are DNA elements that increase transcription of a target gene but can sit thousands to hundreds of thousands of base pairs away, upstream, downstream, or inside introns, and act independently of orientation. They work by physically looping through 3D space to contact the promoter, brought together by architectural proteins (notably CTCF and the cohesin complex). Silencers do the opposite of enhancers — they reduce transcription when bound by repressive factors — and insulators (often CTCF-bound sites) block enhancers from acting on the wrong promoter by restricting which loops can form.

This looping is organized into topologically associating domains (TADs) — megabase-scale genomic neighborhoods, detectable by chromosome conformation capture methods (Hi-C), within which DNA sequences contact each other far more often than they contact sequences in neighboring domains. TAD boundaries, frequently marked by CTCF binding sites, act as soft insulation: an enhancer inside a TAD mostly only reaches promoters inside the same TAD. This is why a genetic variant hundreds of kilobases from a gene can still be its regulator — "nearest gene" is a poor heuristic — and why disrupting a TAD boundary (e.g., by a structural variant) can cause disease by letting an enhancer activate the wrong gene, even with no coding mutation at all.

Failure mode: annotating a variant as "non-coding, likely benign" purely because it falls outside any exon. Non-coding does not mean non-functional; a huge share of disease-associated variants from genome-wide association studies sit in enhancers.

2.6 Chromatin, histone modifications, DNA methylation

DNA does not float free in the nucleus; it wraps around protein spools called histones (eight-protein complexes) to form nucleosomes, and the whole DNA-protein assembly is chromatin. Chromatin exists in a spectrum from open, accessible euchromatin (actively transcribed or transcribable) to tightly packed heterochromatin (silenced). Chemical tags added to histone tails — histone modifications — encode regulatory state without changing the DNA sequence itself, which is why this layer is called epigenetic (literally "on top of" genetics).

DNA methylation adds a methyl group directly to cytosine bases, almost always at CpG dinucleotides (a C followed by a G) in mammals, forming 5-methylcytosine. Promoters with dense clusters of CpGs ("CpG islands") are typically unmethylated when a gene is active and methylated when it is silenced — this is the mechanism behind genomic imprinting and much of tumor-suppressor gene silencing in cancer. Methylation is measured by bisulfite sequencing (which chemically converts unmethylated C to U, leaving methylated C unchanged) and is the molecular layer behind the "methylation array" and "WGBS" data types covered in Module 12 (Epigenomics).

2.7 RNA processing: capping, splicing, alternative splicing, polyadenylation, isoforms

A newly made Pol II transcript (pre-mRNA) is not yet usable; it undergoes three coupled modifications before export to the cytoplasm. First, a modified guanine nucleotide (the 5' cap, 7-methylguanosine) is added to the 5' end, protecting the RNA from degradation and marking it for translation. Second, splicing removes introns (non-coding intervening sequences) and joins exons (the retained, coding or UTR-containing sequences) via the spliceosome, a large RNA-protein machine that recognizes short consensus sequences at splice sites (GT at the intron's start, AG at its end — the "GT-AG rule"). Third, the 3' end is cleaved and a poly(A) tail (a run of adenine nucleotides, typically 50-250 in mature cytoplasmic mRNA) is added, directed by a polyadenylation signal sequence (commonly AAUAAA).

Because exons can be included or skipped combinatorially, a single gene can produce multiple mature mRNA isoforms — alternative splicing. Forms include exon skipping, intron retention, alternative 5'/3' splice sites, and alternative first or last exons (often paired with alternative promoters or polyadenylation sites). The human genome has roughly 20,000 protein-coding genes but is estimated to produce several times that many distinct transcript isoforms, which is why "gene-level" and "transcript-level" quantification in RNA-seq (Module 7) can give materially different biological conclusions from the same data.

Failure mode: quantifying expression at the gene level when the biology of interest is isoform-specific (e.g., a splicing-disrupting variant); gene-level counts can look completely unchanged while isoform usage flips entirely.

2.8 Transcriptome composition and why library prep must choose

"The transcriptome" is not one kind of molecule — it is a mixture dominated, by mass, by ribosomal RNA.

Because rRNA dominates by mass, sequencing total RNA without intervention mostly sequences rRNA. Library preparation protocols must deliberately select for the RNA type of interest:

Failure mode: choosing poly(A) selection and then asking why a known lncRNA or circRNA of interest is absent from the data — it was removed by the protocol, not biologically absent from the cell.

2.9 Translation, the genetic code, codon usage

Translation reads mRNA in non-overlapping triplets called codons, each specifying one amino acid or a stop signal, via the genetic code. The code is redundant (degenerate): 61 codons specify only 20 amino acids, so most amino acids have multiple synonymous codons (e.g., leucine has six). It is read by the ribosome, with transfer RNAs (tRNAs) each carrying a specific amino acid and an anticodon that base-pairs with the mRNA codon. Translation starts at an AUG codon (methionine) within a Kozak consensus context and ends at one of three stop codons (UAA, UAG, UGA).

The code is "near-universal" — a handful of organisms and organelles (notably human mitochondria) reassign a few codons, which matters if you ever translate mitochondrial sequence with the standard nuclear code table and get a frameshift of meaning rather than of sequence.

Organisms do not use synonymous codons equally; codon usage bias reflects tRNA abundance and translation efficiency, measured as the relative synonymous codon usage (RSCU):

$$ RSCU_{i} = \frac{x_i}{\frac{1}{n}\sum_{j=1}^{n} x_j} $$

Here $x_i$ is the observed count of codon $i$, $n$ is the number of synonymous codons for that amino acid, and the denominator is the count expected if all synonyms were used equally; an RSCU above 1 means that codon is preferred relative to chance. This matters practically for heterologous protein expression (codon-optimizing a human gene for expression in E. coli) and for interpreting silent (synonymous) mutations, which can still affect translation speed and protein folding despite not changing the amino acid sequence.

# Translating a CDS with Biopython, respecting the standard genetic code
from Bio.Seq import Seq

cds = Seq("ATGGCCATTGTAATGGGCCGCTGAAAG")
protein = cds.translate(to_stop=True)
print(protein)   # -> Seq('MAIVMGR')

2.10 Protein structure levels, folding, PTMs, domains and families

A protein's primary structure is its linear amino acid sequence, directly determined by the mRNA's codons. Secondary structure describes local, repeating backbone shapes stabilized by hydrogen bonds — mainly the alpha helix and the beta sheet. Tertiary structure is the full three-dimensional fold of a single polypeptide chain, driven largely by hydrophobic amino acid side chains packing into a buried core away from water. Quaternary structure describes the arrangement of multiple folded chains (subunits) into a functional complex, such as hemoglobin's four subunits. Folding happens spontaneously for many proteins following Anfinsen's principle (sequence determines structure under physiological conditions), though many proteins require helper proteins called chaperones to fold correctly in the cell, and misfolding underlies diseases such as the prion disorders mentioned in Section 2.3.

After translation, proteins are frequently chemically modified — post-translational modifications (PTMs) — which expand functional diversity far beyond what the genome alone encodes.

Proteins are also organized into domains — compact, independently folding sequence/structure units that often carry a specific function (a DNA-binding domain, a kinase domain) — and recurring domains across many proteins define protein families (e.g., the kinase family, the immunoglobulin family), cataloged in databases such as Pfam and InterPro. Domain composition is why homology (shared evolutionary origin) can often be inferred from sequence alone, and it is the basis of the functional annotation and protein-family methods covered in Module 10 (Structural Bioinformatics) and Module 11 (Functional Genomics).

Failure mode: assuming sequence identity implies identical function. Two proteins can share a domain and overall fold yet differ in specificity because of a handful of residues in a binding pocket — structure and sequence similarity are evidence for function, not proof of it.

2.9 Gene structure in coordinates: exon, intron, UTR, CDS, strand

A gene is not a single span of DNA with one meaning end to end. When you look at a gene model you are looking at a stack of nested intervals, each with a distinct biological role:

A transcript's exons, read in order, give you the mature RNA. The CDS is a sub-interval of that concatenated exon sequence, not a separate genomic region — this is why "exon 1 is non-coding" is a completely normal statement for a gene whose start codon sits in exon 2.

Strand records which of the two DNA strands is read as the coding (sense) strand for that gene. Genes on the + strand are transcribed left-to-right in increasing coordinate order; genes on the - strand are transcribed right-to-left, so their first exon in transcription order has the highest genomic coordinate. Every tool that reports "upstream" or "downstream," "5' end" or "3' end," means it relative to strand, not relative to the chromosome's coordinate axis. Forgetting this is one of the most common silent bugs in genomics scripts: computing "distance to TSS (transcription start site)" by simple subtraction without checking strand gives the wrong sign for half of all genes.

0-based vs 1-based coordinates — the trap that breaks pipelines

Different file formats count genomic positions differently, and this is not a stylistic footnote — it silently shifts every interval by one base if you mix formats without converting.

• 1-based, closed (inclusive) interval: position 1 is the first base of the chromosome; an interval [start, end] includes both endpoints. Used by: GFF3, GTF, VCF, SAM, most genome browsers' on-screen coordinates.
• 0-based, half-open interval: position 0 is the first base; an interval [start, end) includes start but excludes end. Used by: BED, BAM (internally), most Python/C-style array logic.

Worked example. Suppose you want to describe the third and fourth bases of a chromosome, i.e., bases at 1-based positions 3 and 4.

Both rows describe the exact same two bases. The conversion rule is: BED_start = GFF_start - 1, BED_end = GFF_end (end is unchanged because closed-inclusive end in 1-based equals open-exclusive end in 0-based when you've already shifted the start down by one). Get this backwards and every feature in your BED file is shifted by one base relative to GFF-derived coordinates — a single-base error that is invisible on a genome browser zoomed out, catastrophic for variant-to-exon overlap calls at splice sites, and a classic source of irreproducible results when two tools (one BED-based like bedtools, one GFF-based like a GTF parser) are combined without conversion.

# bedtools intersect expects BED-style (0-based, half-open) input
# If you build a BED record "by hand" from a 1-based GFF coordinate, subtract 1 from start only:
# GFF: chr1  5  10   (1-based, inclusive: bases 5-10)
echo -e "chr1\t4\t10" > region.bed   # correct BED equivalent: 0-based start, end unchanged

A practical rule: never hand-edit coordinates between formats. Use a library that knows the convention (pybedtools, pyranges, GenomicRanges in Bioconductor) and let it do the arithmetic.

2.10 Genome assemblies and reference versions

A coordinate by itself means nothing. "Chromosome 7, position 140,453,136" is a different base depending on which reference assembly you mean, because successive assemblies correct errors, fill gaps, and renumber regions.

"hg19" and "GRCh37" are, for almost all practical purposes, the same coordinate system (there are minor chromosome-naming and a few sequence-patch differences, chr1 vs 1 being the most visible). Mixing hg19-called variants with a GRCh38 annotation file without converting is a guaranteed way to silently place variants in the wrong gene, because the same numeric position refers to different DNA.

Liftover is the process of converting coordinates from one assembly to another using a precomputed alignment ("chain file") between the two assemblies.

# UCSC liftOver: convert hg19 BED intervals to hg38
liftOver input_hg19.bed hg19ToHg38.over.chain.gz output_hg38.bed unmapped.bed
# unmapped.bed lists intervals that could not be confidently mapped (duplicated, deleted, or rearranged regions)

# CrossMap is the equivalent tool for VCF, BAM, GFF, and BigWig, not just BED
CrossMap.py vcf hg19ToHg38.over.chain.gz input_hg19.vcf GRCh38.fa output_hg38.vcf

Liftover is not always exact: regions that were rearranged, duplicated, or deleted between assembly versions cannot be mapped one-to-one, and a fraction of intervals (typically well under 1% genome-wide, but locally much higher in repetitive or structurally variable regions) will fail or map ambiguously. For clinical or publication-grade work, re-calling variants directly against the target assembly is safer than lifting over calls made on another assembly. Always record the assembly version in every file name, header line, and figure legend — a BED file or VCF with no assembly tag is not reusable six months later, including by the person who made it.

2.11 Identifiers: symbols, Ensembl, RefSeq, UniProt, Entrez — and why this is a minefield

The same biological entity accumulates multiple identifiers from different databases, and these identifiers are not interchangeable or always one-to-one.

Why symbols are dangerous. A gene symbol is a human-readable label, not a database key, and it fails in three specific ways:

1. One symbol, multiple meanings across time or species. Symbols get reassigned by HGNC as understanding of gene families improves; an old paper's "KIAA0101" is now "PCLAF." Mouse and human orthologs often share a symbol differing only in case (Tp53 vs TP53), which breaks naive string joins.
2. One gene, multiple symbols (aliases). TP53 is also referenced historically as P53. A gene list built by scraping text will contain a mix of current symbols and aliases that a simple lookup table will not reconcile unless it explicitly includes alias mapping.
3. The Excel gene-symbol disaster. Spreadsheet software auto-converts certain gene symbols into dates or scientific notation on file open: SEPT1 (Septin 1) becomes 1-Sep, MARCH1 becomes 1-Mar, and DEC1 can become a numeric date serial. A peer-reviewed survey of supplementary gene lists found that roughly one in five genomics papers with Excel supplementary tables had this corruption somewhere in the file. The fix is procedural, not clever: never open a gene list in Excel with default settings; import gene-symbol columns explicitly as text, or avoid spreadsheet software for data interchange entirely and use TSV/CSV read programmatically.

ID mapping pitfalls in practice. Gene-to-transcript-to-protein is one-to-many in both directions once you account for alternative splicing (Module 2.6) — one gene has several transcripts, and some transcripts' CDS can, rarely, map to the same protein sequence. Cross-database mapping tables (Ensembl BioMart, NCBI's gene2ensembl, UniProt's ID mapping service) are the correct tool, but they still drop or multiply rows silently if you don't check cardinality after every join.

import mygene

mg = mygene.MyGeneInfo()
result = mg.querymany(
    ["TP53", "SEPT1", "BRCA1"],
    scopes="symbol",
    fields=["ensembl.gene", "entrezgene", "uniprot.Swiss-Prot"],
    species="human",
)
# Inspect result for 'notfound': True entries and for queries that
# returned more than one hit (ambiguous symbol) before trusting any join.
for r in result:
    print(r.get("query"), r.get("ensembl", {}).get("gene"), r.get("notfound"))

Always version-pin Ensembl IDs (ENSG00000141510.18, not just ENSG00000141510) in any pipeline you expect to reproduce, because the gene model behind that ID — its exon structure, its canonical transcript — can change between Ensembl releases even though the stable ID stays the same.

2.12 Genetic variation: the taxonomy

The size thresholds separating these categories (50 bp for indel/SV boundary, 1 kb-ish for CNV) are conventions used by variant callers, not biological laws — different tools and consortia draw the line slightly differently, so always check a given caller's or database's definition before comparing counts across studies.

2.13 Allele frequency, genotype, zygosity, haplotype, linkage disequilibrium, Hardy-Weinberg

An allele is one version of a sequence at a given locus. A genotype is the pair of alleles an individual carries at that locus (diploid organisms have two). Zygosity describes that pair: homozygous (two identical alleles), heterozygous (two different alleles), or hemizygous (only one copy present, as for X-linked genes in males). A haplotype is a set of alleles at multiple linked loci inherited together on the same physical chromosome, because recombination between very close loci is rare.

Allele frequency is simply the proportion of all alleles in a population that are a given variant. Linkage disequilibrium (LD) measures whether two alleles at different loci co-occur more or less often than chance would predict given their individual frequencies, commonly summarized as $r^2$ (ranges 0 to 1, where 1 means the two loci are perfectly predictive of each other) or $D'$.

Hardy-Weinberg equilibrium is the expected genotype distribution in a large, randomly mating population with no selection, mutation, or migration acting on a locus:

$$p^2 + 2pq + q^2 = 1$$

Here $p$ is the frequency of one allele and $q = 1-p$ is the frequency of the other; $p^2$ is the expected fraction of homozygotes for the first allele, $q^2$ the expected fraction of homozygotes for the second, and $2pq$ the expected fraction of heterozygotes — the factor of 2 appears because a heterozygote can arise from either parent contributing either allele. Large deviations from these proportions in observed genotype counts flag genotyping error, population stratification, or real selection at that locus, which is why a Hardy-Weinberg equilibrium test (chi-square goodness-of-fit) is a standard genotyping quality-control step before any association analysis (Module 10, Statistical Genetics).

from scipy.stats import chisquare

# Observed genotype counts: AA, Aa, aa
obs = [450, 420, 130]
n = sum(obs)
p = (2*obs[0] + obs[1]) / (2*n)   # allele frequency of A
q = 1 - p
exp = [p**2 * n, 2*p*q * n, q**2 * n]
stat, pval = chisquare(obs, exp, ddof=1)   # ddof=1: one parameter (p) estimated from the data
print(p, exp, pval)

2.14 Germline vs somatic, mosaicism, clonality

A germline variant is present in the fertilized egg and therefore in every cell of the body, and is heritable. A somatic variant arises after fertilization, in one cell lineage only, and is not transmitted to offspring; most cancer-driving mutations are somatic. Mosaicism is the presence of two or more genetically distinct cell populations in one individual arising from a post-zygotic mutation — it can be confined to a tissue (segmental mosaicism) or scattered at low level across blood, as in clonal hematopoiesis, where a single hematopoietic stem cell's mutation expands to detectably contribute to the blood cell pool with age. Clonality describes whether a population of cells (a tumor, an expanded immune cell population) descends from one ancestral cell; it is read out computationally from shared mutations or variant allele frequency (VAF), the fraction of sequencing reads supporting a variant, which should cluster near 50% for a heterozygous germline variant but can take any value for somatic variants depending on tumor purity and subclone size. Distinguishing germline from somatic is the single most consequential analytical decision in cancer genomics (Module 15), because it requires a matched normal (non-tumor) sample as comparison — calling variants from tumor tissue alone cannot separate the two.

2.15 Population structure and ancestry

Human populations differ systematically in allele frequencies because of shared ancestry, migration history, and genetic drift, not because of any single deterministic trait. This population structure matters computationally for two reasons: it confounds association studies (a frequency difference between cases and controls can reflect ancestry rather than disease biology) and it determines which reference panel is appropriate for interpreting a given individual's variants.

Principal component analysis (PCA) on genome-wide genotype data is the standard way to visualize and correct for structure: the first few principal components typically separate major continental ancestry groups, and including them as covariates in association models is routine. $F_{ST}$ quantifies differentiation between two populations as the proportion of total genetic variance attributable to allele-frequency differences between them, ranging from 0 (identical frequencies) to 1 (complete separation). Public reference panels — the 1000 Genomes Project, the Human Genome Diversity Project, and allele-frequency databases like gnomAD — provide population-stratified frequencies used to flag whether a variant is common in some populations and rare or absent in others, which is essential context: a variant's "rarity," and therefore its presumed pathogenicity, is only meaningful relative to a matched ancestry background, and tools or clinical pipelines trained predominantly on European-ancestry reference data systematically misestimate frequency and risk for underrepresented populations.

2.18 What a gene's "function" means, and the ontologies that encode it

Saying a gene "does X" is shorthand for several different claims that get blurred in casual speech. A gene can have a molecular function (what the protein does biochemically — "ATP binding", "DNA-binding transcription factor activity"), a role in a biological process (what happens at the cell or organism level — "apoptotic process", "glycolysis"), and a cellular component localization (where it acts — "mitochondrial inner membrane", "nucleus"). These three axes are deliberately kept separate because a protein can have one molecular function but participate in many processes depending on context (cell type, developmental stage, disease state), and the same process can be carried out by different proteins in different species.

The Gene Ontology (GO) formalizes this as three structured, controlled vocabularies (Molecular Function, Biological Process, Cellular Component), each a directed acyclic graph (DAG) — a tree-like structure where a term can have more than one parent. "Mitochondrial inner membrane" is a child of both "mitochondrion" and "membrane". Terms are connected by relationships, mainly is_a (strict subtype) and part_of (component relationship, not subtype — a ribosome is part_of a cell but is not is_a a cell). Every gene is annotated to one or more terms, each annotation carrying an evidence code that tells you how reliable the claim is: EXP (direct experimental evidence), IDA (inferred from direct assay), down to IEA (inferred from electronic annotation — a computational guess, never manually checked). Treating an IEA annotation with the same confidence as an IDA one is a common and avoidable error; always check the evidence code before trusting an enrichment result.

GO enrichment analysis asks: given a list of genes (say, the 300 genes significantly upregulated in an experiment, from Module 7 differential expression), are any GO terms over-represented relative to what you'd expect by chance? The standard test is the hypergeometric test (equivalent to Fisher's exact test on a 2x2 table):

$$P(X \ge k) = \sum_{i=k}^{\min(K,n)} \frac{\binom{K}{i}\binom{N-K}{n-i}}{\binom{N}{n}}$$

Here $N$ is the total number of genes in the background (usually all genes tested, not the whole genome), $K$ is the number of genes annotated to the term of interest, $n$ is the size of your gene list, and $k$ is how many of your genes fall in that term. The formula counts, among all possible ways to draw $n$ genes from $N$, how many draws would give you $k$ or more "hits" in a category of size $K$ — if that number is small, your observed overlap is unlikely to be chance. Because you test thousands of GO terms simultaneously, you must correct for multiple testing (Benjamini-Hochberg false discovery rate; see Module 7) — an uncorrected p-value of 0.01 across 15,000 terms will produce hundreds of false positives.

Pathway databases add the piece GO deliberately omits: order and interaction. GO tells you a gene is "involved in glycolysis"; a pathway resource tells you it catalyzes step 3, consuming this substrate and producing that product, feeding into this other pathway.

A practical distinction: GO and MSigDB give you gene sets (unordered bags of genes with a shared label) suitable for set-overlap tests; KEGG and Reactome give you pathway topology (who talks to whom) suitable for tools that use network structure (e.g., SPIA, Reactome's own pathway browser) to ask not just "are these genes involved" but "is the flow through the pathway perturbed". A minimal enrichment call in R, using a popular wrapper:

library(clusterProfiler)
library(org.Hs.eg.db)

# gene_list: character vector of Entrez IDs, e.g. c("7157","1017","672")
ego <- enrichGO(gene          = gene_list,
                universe      = background_genes,   # all genes tested, not whole genome
                OrgDb         = org.Hs.eg.db,
                ont           = "BP",                # Biological Process
                pAdjustMethod = "BH",
                pvalueCutoff  = 0.05,
                qvalueCutoff  = 0.2)
head(as.data.frame(ego))
# Columns: ID, Description, GeneRatio, BgRatio, pvalue, p.adjust, qvalue, geneID, Count

The universe argument is the most commonly skipped and most consequential argument in this call: if you leave it as the default (all annotated genes genome-wide) instead of the actual set of genes your experiment could have detected (e.g., only expressed genes passing a filter), you will get systematically inflated significance, because your background no longer matches your assay.

2.19 Model organisms and orthology

Most of what we know about gene function comes not from humans but from a small set of model organisms chosen for practical reasons — fast generation time, cheap husbandry, genetic tractability, and decades of accumulated community knowledge. The function annotations in GO, the pathway diagrams in KEGG, and the phenotypes linked to genes in disease databases are heavily built from this cross-species work, so understanding how function is transferred between species is essential to reading any annotation correctly.

Orthologs are genes in different species descended from a single gene in their most recent common ancestor via speciation; they typically retain similar function. Paralogs arise from a gene duplication event, within one lineage; they may diverge in function. This distinction matters enormously for inference: if human gene X has no direct functional data but its mouse ortholog has been knocked out and characterized, you can reasonably infer function by orthology. If you instead find a human paralog of X that diverged before a key duplication, inferring the same function is unjustified — paralogs are notorious for having drifted into different roles (classic case: the globin gene family — hemoglobin subunits and myoglobin are paralogs with related but distinct jobs).

Orthology is not always one-to-one. One-to-one orthology (one gene, one counterpart) is the clean case. One-to-many and many-to-many orthology occur when lineage-specific duplications happened after the species split — common in gene families like olfactory receptors or immune gene clusters, where the mouse may have five paralogous copies of a gene that exists as one copy in human. Tools that compute orthology (Ensembl Compara, OrthoDB, PANTHER, the DIOPT meta-tool that aggregates several methods) report a confidence level, and you should treat a "many-to-many, low confidence" call as a hypothesis, not a fact. A frequent pipeline mistake is naively mapping gene symbols across species 1:1 (Trp53 in mouse to TP53 in human by string similarity) without consulting an orthology database — this silently fails for the thousands of genes with non-obvious name correspondence or genuine multi-copy relationships.

2.20 Cells, tissues, development, and differentiation

A genome is a single, mostly fixed set of instructions, but a human body contains several hundred distinct cell types, every one carrying (with minor exceptions — see somatic mosaicism in the earlier section on variation) the identical DNA sequence. The difference between a neuron and a hepatocyte is not genomic, it is in which genes are switched on and read, mediated by the chromatin state, transcription factor networks, and epigenetic marks covered in sections 2.5-2.6. This is the central fact that single-cell and spatial genomics (Modules 12-13) are built to measure: not "what is in the genome" but "what is each individual cell actually doing right now, and where is it".

Development is the process by which a single fertilized cell (the zygote) divides and diversifies into all these cell types, following a program of sequential gene-expression decisions. Differentiation is the narrowing of a cell's potential — a stem cell (capable of both self-renewal and producing differentiated progeny) becomes progressively more specialized. Potency is described on a gradient:

A tissue is an organized assembly of cell types (plus extracellular matrix, vasculature, immune infiltrate) cooperating in a function; an organ is one or more tissues arranged into a larger functional unit. The same cell type can look and behave differently depending on tissue context (a fibroblast in skin vs. lung), which is why cell type is increasingly defined operationally by its molecular state — a transcriptomic or epigenomic signature — rather than by morphology alone.

This reframing is exactly what single-cell RNA-seq (scRNA-seq, Module 12) and spatial transcriptomics (Module 13) exploit. Bulk RNA-seq (Module 7) measures average expression across millions of cells in a sample and cannot tell you whether a signal comes from a few cells expressing a gene strongly or many cells expressing it weakly — this is cellular heterogeneity, and it is invisible to bulk assays by construction. scRNA-seq profiles thousands of individual cells, clusters them by expression similarity, and assigns each cluster a cell-type label using marker genes (genes known from prior literature to be specific to one type, e.g., PTPRC/CD45 for all leukocytes, EPCAM for epithelial cells). Spatial methods go one step further, retaining each cell's (or small region's) physical coordinates in the tissue, so you can ask not just "which cell types are present" but "which cell types are next to each other" — directly relevant to questions like whether immune cells are infiltrating a tumor or being excluded from it. The biology in this section is the prerequisite for interpreting every UMAP plot and every spatial cluster map you will encounter later in the course: a cluster is only meaningful insofar as it corresponds to a real, biologically coherent cell state, and that correspondence has to be argued with marker genes and known developmental biology, not assumed from the plot alone.

2.21 Cancer as a disease of the genome

Cancer is fundamentally a disease in which cells accumulate genomic and epigenomic changes that let them escape the normal controls on proliferation, survival, and tissue boundaries. It is not one disease but a category defined by this shared mechanism, which is why "the genomics of lung cancer" and "the genomics of leukemia" can look almost unrelated at the level of which specific genes are altered, while sharing the same underlying logic of cause.

The key conceptual vocabulary:

• Oncogene: a gene which, when activated (by mutation, amplification, or translocation), drives proliferation. Typically requires only one mutant allele to act (dominant). Example: KRAS point mutations (G12D, G12C) that lock the protein in an active, growth-signaling state.
• Tumor suppressor gene: a gene whose normal job is to restrain proliferation or trigger repair/death of damaged cells; cancer typically needs both copies disabled (recessive at the cell level) — Knudson's "two-hit hypothesis". Example: TP53, mutated or deleted in roughly half of all human cancers; BRCA1/BRCA2, whose germline loss predisposes to breast and ovarian cancer.
• Driver vs. passenger mutation: a driver mutation confers a selective growth advantage to the cell carrying it and is causally implicated in cancer progression; a passenger mutation is along for the ride, acquired incidentally during the many cell divisions of tumor growth, with no functional consequence. A typical solid tumor carries tens to thousands of somatic mutations (see Module 10/15 on somatic variant calling), and distinguishing the handful of drivers from the passenger background is one of cancer genomics' central statistical problems, usually tackled by recurrence across patients (is this exact mutation, or this gene, hit far more often than the background mutation rate would predict) rather than by inspecting a single tumor.
• Clonal evolution: a tumor is not a static genome but a population under selection. It starts from one cell, accumulates mutations, and the fittest subclones expand, producing intratumor heterogeneity — a single tumor can contain multiple genetically distinct subclones (connects directly to the clonality concept introduced earlier in this module for somatic mosaicism). Metastases and treatment resistance are frequently traced to the expansion of a pre-existing minor subclone under the selective pressure of therapy.
• Tumor mutational burden (TMB): the number of somatic mutations per megabase of tumor genome. High-TMB tumors tend to present more neoantigens (mutant peptides that can be recognized as foreign) to the immune system and are more likely to respond to immune checkpoint inhibitor therapy — a direct, clinically actionable link between genomic data and treatment choice.

Structural events matter as much as point mutations. Copy number amplification of an oncogene (e.g., ERBB2/HER2 amplification in breast cancer, which is the direct basis for trastuzumab therapy eligibility) increases dosage of an already-active gene. Gene fusions from chromosomal translocation can create a novel, constitutively active protein — the BCR-ABL1 fusion from the Philadelphia chromosome translocation t(9;22) is the textbook case, and it is also a direct drug target (imatinib). Loss of heterozygosity (LOH) — loss of the wild-type copy of a gene whose other copy was already mutated — is the common mechanistic route to Knudson's second hit. Each of these event types requires different computational detection strategies covered in Modules 10 and 15: point mutations need deep, paired tumor-normal sequencing and careful somatic-vs-germline filtering; copy number needs read-depth and allele-balance analysis; fusions need split-read and discordant-read detection or RNA-seq-based callers.

2.22 The immune system in one page

Immunology shows up in essentially every modern dataset — tumor biopsies contain immune infiltrate, single-cell atlases are dominated by immune cell diversity, vaccine and infectious-disease studies are immunology directly — so a minimal working vocabulary is necessary even outside immunology-specific work.

The immune system has two broad arms. Innate immunity is fast, non-specific, and present from birth: physical barriers, phagocytic cells (macrophages, neutrophils) that engulf pathogens, and pattern-recognition receptors that detect generic molecular signatures of infection (bacterial cell wall components, viral RNA). It acts within minutes to hours and does not improve with repeated exposure. Adaptive immunity is slow to initiate (days) but highly specific and carries memory — repeated exposure to the same pathogen produces a faster, stronger response. Its two main cell types are T cells, which mature in the thymus and either kill infected cells directly (cytotoxic/CD8+ T cells) or coordinate other immune cells (helper/CD4+ T cells), and B cells, which mature in the bone marrow and, upon activation, differentiate into plasma cells that secrete antibodies (proteins that bind a pathogen and tag it for destruction or neutralize it directly).

The specificity of adaptive immunity rests on two independent diversity-generating mechanisms that are directly relevant to sequencing data:

• MHC/HLA: the Major Histocompatibility Complex (called HLA, Human Leukocyte Antigen, in humans) is a set of genes encoding proteins that sit on the cell surface and display short peptide fragments from inside the cell (or from engulfed material) to T cells. HLA genes are the most polymorphic loci in the human genome — thousands of alleles exist at genes like HLA-A, HLA-B, HLA-DRB1 — because different alleles present different peptides well, and population-level diversity hedges against any single pathogen evading detection entirely. This is why organ and bone marrow transplant matching requires HLA typing, why some drug hypersensitivity reactions are HLA-allele-specific (e.g., abacavir and HLA-B*57:01), and why standard short-read variant callers often misassign reads in the HLA region — its sequence diversity and structural complexity routinely break assumptions about a single linear reference, requiring specialized HLA-typing tools (Module 10) rather than generic variant callers.
• TCR/BCR diversity: T cell receptors (TCR) and B cell receptors (BCR, the membrane-bound precursor of an antibody) achieve enormous diversity (theoretically beyond $10^{15}$ distinct sequences) through V(D)J recombination — each receptor gene is assembled during cell development by randomly selecting and joining one Variable, one (for some chains) Diversity, and one Joining gene segment from a repertoire of options, with additional random nucleotide insertion/deletion at the joints. Every mature T or B cell carries one fixed receptor sequence, and a population of cells sharing an identical recombined sequence is a clone, descended from one ancestor cell that first encountered its matching antigen. Immune repertoire sequencing (a specialized library prep and analysis track touched on in Module 12) reads out exactly this: which V(D)J combinations are present and at what frequency, used to track immune responses to infection, vaccination, and tumors, and to detect minimal residual disease in lymphoid cancers by tracking a known malignant clone's exact receptor sequence.

2.23 Biological object to data representation — the master map

2.24 Common pitfalls and how to avoid them

2.25 Exercises

1. (Warm-up) Define, in one sentence each, the difference between an ortholog and a paralog, and give one worked example of each from this section. Deliverable: three sentences.
2. (Warm-up) Using the master table in section 2.23, state which file format you would receive from a single-cell RNA-seq experiment and which from a Hi-C (3D genome contact) experiment, and name the module where each is covered in depth. Deliverable: a two-row table.
3. (Core) You run enrichGO on 250 upregulated genes from an RNA-seq experiment that tested 14,000 expressed genes, using the default genome-wide universe of ~20,000 genes instead of your 14,000 tested genes. Explain, in mechanistic terms referencing the hypergeometric formula in section 2.18, why this inflates significance. Deliverable: a short paragraph identifying which term in the formula is wrong and in which direction.
4. (Core) A tumor sample shows a TP53 mutation in one allele and a copy-number loss spanning the TP53 locus on the other chromosome. Using the vocabulary from section 2.21, name this combination of events and explain why it is sufficient to inactivate a tumor suppressor gene when a single point mutation alone would not be. Deliverable: a short paragraph naming both the "two-hit" concept and the LOH mechanism.
5. (Core) Explain why a cell biologist studying HLA genotyping would not simply run the same germline variant-calling pipeline used for the rest of the genome. Deliverable: 3-4 sentences referencing HLA polymorphism and reference-bias.
6. (Stretch) Design, in prose (no code required), an experiment distinguishing whether an observed difference in average expression of gene X between two tissue samples is due to (a) every cell changing its expression level uniformly, or (b) a change in the proportion of a cell type that already expresses gene X highly. Deliverable: a short experimental plan naming which assay type (bulk vs. single-cell) resolves the ambiguity and why.
7. (Stretch) A TCR repertoire sequencing experiment on a patient's tumor-infiltrating lymphocytes finds one CDR3 sequence present at 40% of all reads, while the rest of the repertoire is highly diverse. Propose a biological interpretation consistent with section 2.22's vocabulary, and state what additional data (e.g., matched blood sample, antigen specificity assay) would support or refute it. Deliverable: one paragraph.

Solutions / hints

1. Orthologs arise by speciation from a shared ancestral gene (example: human TP53 and mouse Trp53); paralogs arise by gene duplication within a lineage and can diverge in function (example: the hemoglobin and myoglobin globin-family paralogs).
2. Single-cell RNA-seq: cell-by-gene sparse matrix, typically .h5ad (AnnData) or Seurat .rds, Module 12. Hi-C: pairwise contact frequencies, .hic or cooler .cool/.mcool, Module 9.
3. The formula's $N$ (background size) is wrongly set to ~20,000 instead of 14,000, and $K$ (genes annotated to a term) is computed over that same wrong, larger background. Because many genes in the extra ~6,000 are not expressed in this tissue/assay at all, the apparent "non-hit" background pool is artificially inflated with genes that could never have appeared in your list, which shrinks the computed p-value for any observed overlap — making enrichment look stronger (more significant) than it really is relative to the genes you actually could have detected.
4. This is Knudson's two-hit scenario: the point mutation inactivates one TP53 allele, and the copy-number loss removes the other copy entirely — a form of loss of heterozygosity (LOH). Because tumor suppressors are typically recessive at the cellular level, one functional copy is enough to maintain normal restraint on proliferation; only when both copies are disabled (by any combination of mutation, deletion, or epigenetic silencing) is the brake fully released.
5. HLA genes are the most polymorphic loci in the human genome, with thousands of alleles; a standard pipeline aligns reads to one reference haplotype per gene, so reads from divergent HLA alleles may fail to align, align with excess mismatches, or map to the wrong HLA gene entirely, biasing both read depth and genotype calls. Dedicated HLA callers use curated allele databases and specialized alignment/EM-based allele-assignment algorithms instead of a single linear reference.
6. Bulk RNA-seq cannot distinguish (a) from (b) because it only reports the population average. Single-cell RNA-seq resolves it directly: cluster cells by type, and check whether the per-cell expression level of gene X within the matching cell type differs between samples (supports a) versus whether the proportion of high-expressing cells differs while per-cell level is stable (supports b).
7. A CDR3 sequence at 40% of reads indicates a dominant T cell clone, consistent with a strong, antigen-driven expansion of one T cell lineage — plausibly an anti-tumor response if the antigen is tumor-derived. To support or refute this, compare against a matched peripheral blood sample (is the clone tumor-restricted or systemic?) and, if possible, test the clone's TCR for specificity against candidate tumor antigens (e.g., via a tetramer assay or functional antigen-stimulation assay).

2.26 Key takeaways

• A gene's "function" is really three separable claims (molecular function, biological process, cellular location); GO keeps them as three ontologies for this reason.
• GO terms form a DAG, not a tree — a term can have multiple parents via is_a and part_of relationships.
• Always check the evidence code behind a GO annotation; IEA is a computational guess, not a validated fact.
• Enrichment analysis background (universe) must match the genes your assay could actually detect, not the whole genome.
• KEGG and Reactome encode pathway wiring; GO and MSigDB encode unordered gene sets — pick the resource that matches the question.
• Orthologs (shared ancestor via speciation) usually conserve function; paralogs (shared ancestor via duplication) often do not — never assume 1:1 symbol matching across species.
• Model organisms underlie most functional annotation in human gene databases; the inference chain from model organism to human is only as strong as the orthology call.
• Differentiation is a change in gene expression program, not genome sequence; single-cell and spatial assays exist specifically to resolve this cell-by-cell.
• Cancer genomics separates driver mutations (causal, recurrent) from passenger mutations (incidental); tumor mutational burden links directly to immunotherapy response.
• Tumor suppressors typically need both alleles disabled (Knudson two-hit, often via loss of heterozygosity); oncogenes typically need only one activating event.
• HLA diversity and V(D)J recombination are the two independent engines of immune specificity, and both break standard reference-based bioinformatics assumptions.
• Every biological object in this module maps to a specific file format and a specific later module — this map is the skeleton the rest of the course hangs on.

2.27 Further reading

• Gene Ontology Consortium — official site and documentation: geneontology.org; foundational paper: Ashburner et al., "Gene Ontology: tool for the unification of biology".
• KEGG (Kyoto Encyclopedia of Genes and Genomes) — kegg.jp; Kanehisa and Goto, original KEGG database paper.
• Reactome — reactome.org, open-access curated pathway database.
• WikiPathways — wikipathways.org, community-curated pathway database.
• MSigDB (Molecular Signatures Database) — Broad Institute GSEA/MSigDB site; Subramanian et al., "Gene set enrichment analysis" (original GSEA paper).
• clusterProfiler R package documentation (Bioconductor) for GO/KEGG enrichment workflows.
• Ensembl Compara and OrthoDB documentation for orthology/paralogy inference methods.
• DIOPT (DRSC Integrative Ortholog Prediction Tool), Harvard Medical School, for cross-database ortholog consensus calls.
• International Mouse Phenotyping Consortium (IMPC) — mousephenotype.org, systematic mouse knockout phenotyping.
• COSMIC (Catalogue of Somatic Mutations in Cancer) and its Cancer Gene Census — curated list of known cancer driver genes.
• Knudson, "Mutation and cancer: statistical study of retinoblastoma" — original two-hit hypothesis paper.
• IPD-IMGT/HLA Database — official reference for HLA allele nomenclature and sequences.
• AIRR Community (Adaptive Immune Receptor Repertoire Community) — standards body for TCR/BCR repertoire data formats.
• Janeway's Immunobiology (textbook) for a fuller treatment of innate/adaptive immunity than this one-page summary.

Module 3 — Sequences, Alignment, and the Databases

In one paragraph. This module treats biological sequences as what they are to a computer: strings over a small alphabet, stored in specific file formats with specific rules, and compared using dynamic-programming algorithms that are over 50 years old but still underlie every modern aligner. You will learn to read and reason about FASTA, FASTQ, SAM/BAM/CRAM, VCF/BCF, BED/GFF3/GTF, and the common signal-track and single-cell container formats field by field, not by vague description. You will also work two classic alignment algorithms by hand on paper-sized examples, then implement them in Python, so that when a production tool like BWA or minimap2 gives you an unexpected CIGAR string or a MAPQ of zero, you know exactly what happened inside the box.

Prerequisites: Module 1 (orientation: what DNA/RNA/protein are, central dogma), Module 2 (command line, basic Python, working with files and environments). No prior bioinformatics-specific knowledge assumed. You will be able to: - Compute reverse complements, GC content, and k-mer spectra from raw sequence, and explain why repeats break assembly and mapping. - Read and write FASTA and FASTQ correctly, including quality encodings and Illumina header fields. - Decode any SAM FLAG value, parse any CIGAR string by hand, and explain what MAPQ actually measures. - Distinguish INFO from FORMAT fields in a VCF and interpret GT/AD/DP/GQ/PL for a real genotype call. - Explain the coordinate-system difference between BED and GFF3/GTF and convert between them without off-by-one errors. - Derive, conceptually, where a BLOSUM or PAM substitution matrix's numbers come from. - Hand-compute a Needleman-Wunsch and a Smith-Waterman alignment matrix and traceback. - Implement global and local alignment in Python from scratch, and reproduce the same result with Biopython, parasail, or edlib.

Time: 6-8 hours (longer if you work every hand-computation yourself, which you should).

3.1 Sequence data as strings

A DNA sequence is, to every piece of software you will ever run, a string over a finite alphabet. Biology adds meaning on top, but the file on disk is text (or a tightly packed binary encoding of text). Getting comfortable with string-level thinking is what lets you debug pipelines instead of treating them as oracles.

Alphabets. DNA uses {A, C, G, T} (and N for "unknown base"). RNA replaces T with U. Protein uses the 20 standard amino acid letters plus a handful of ambiguity and special codes (X = any amino acid, * or . = stop codon, U = selenocysteine, O = pyrrolysine in a few organisms). DNA sequencing and genotyping also use the IUPAC ambiguity codes, which represent "the true base is one of these two or more possibilities" — essential for representing heterozygous sites or degenerate primers.

These codes appear in primer design (a degenerate primer targeting a polymorphic site), in reference genomes at masked or unresolved positions, and occasionally in variant call output for symbolic alleles. Any parser you write for raw sequence must not silently assume only ACGT — real reference genomes (e.g. GRCh38) contain long stretches of N for unresolved centromeric/telomeric gaps, and real reads can contain N where the sequencer had no confidence in a base call.

Reverse complement. DNA is double-stranded and antiparallel: the strand you did not sequence is the complement of the one you did, read 3′→5′ relative to it, which conventionally we always write 5′→3′, hence "reverse complement." Complementation pairs A↔T and C↔G (Watson-Crick base pairing: two hydrogen bonds for A-T, three for G-C, which is why GC-rich DNA is thermally and mechanically more stable). You need the reverse complement constantly: to check whether a read came from the minus strand, to design a primer that binds the opposite strand, to search a sequence database for both orientations.

COMPLEMENT = str.maketrans("ACGTNacgtn", "TGCANtgcan")

def reverse_complement(seq: str) -> str:
    return seq.translate(COMPLEMENT)[::-1]

reverse_complement("GATTACA")
# 'TGTAATC'

GC content. The fraction of bases that are G or C:

$$\mathrm{GC} = \frac{n_G + n_C}{n_A + n_C + n_G + n_T}$$

where $n_X$ is the count of base $X$. GC content varies by organism (as low as ~20% in some AT-rich parasites, over 70% in some soil bacteria) and even along a single genome ("isochores" in vertebrates). It matters practically because PCR primers, sequencing coverage uniformity, and even assembly difficulty are GC-content-dependent: PCR under- or over-amplifies extreme-GC regions, and Illumina coverage famously dips in both very GC-poor and very GC-rich windows due to amplification bias during library prep.

k-mers. A k-mer is every substring of length $k$ that occurs in a sequence, read with a sliding window of stride 1. A sequence of length $L$ has $L - k + 1$ k-mers (fewer than that many distinct k-mers, since some repeat). K-mers are the currency of modern genome-scale algorithms precisely because fixed-length strings can be hashed, sorted, and counted in near-linear time, whereas full alignment is quadratic. De Bruijn graph assemblers (Module 4) build a graph whose nodes are k-mers; read classifiers like Kraken2 and error-correctors like BFC compare k-mer spectra against reference sets; mash and sourmash estimate genome similarity (distance) from k-mer sketches without ever aligning anything.

def kmers(seq: str, k: int) -> list[str]:
    return [seq[i:i+k] for i in range(len(seq) - k + 1)]

kmers("GATTACA", 3)
# ['GAT', 'ATT', 'TTA', 'TAC', 'ACA']

A k-mer spectrum is the histogram of how many times each distinct k-mer occurs across a dataset. For a correctly-sequenced genome at uniform coverage $c$, most true k-mers occur approximately $c$ times (a Poisson-ish peak), and sequencing-error-induced k-mers occur once or twice — this is exactly the signal used to estimate genome size and heterozygosity (tools like GenomeScope) and to flag likely sequencing errors before assembly.

Complexity and low-complexity sequence. Not all strings are equally "surprising." A run like AAAAAAAAAA or a short tandem repeat like CACACACACA has low sequence complexity: it can be described (compressed) much more briefly than its literal length, because a short pattern repeats. One common formal measure is Shannon entropy over a sliding window, computed from the frequency $p_i$ of each symbol:

$$H = -\sum_i p_i \log_2 p_i$$

where $p_i$ is the fraction of positions in the window occupied by symbol $i$. $H$ is maximal (2 bits for 4-letter DNA) when all bases are equally frequent and locally unpredictable, and drops toward 0 as the window becomes dominated by one or two repeating symbols. Low-complexity sequence matters because standard alignment scoring (Section 3.3) rewards matches — a region of AAAAA will align to almost any other AAAAA-containing sequence with a high score that reflects nothing biological, it is a statistical trap. Tools like dustmasker, RepeatMasker, and the built-in low-complexity filters in BLAST exist specifically to soft-mask (lowercase) or hard-mask (replace with N) these regions before alignment, so that downstream hit-scoring is not dominated by compositional bias.

Repeats, and why they break everything. A repeat is any sequence that occurs more than once in a genome, at near-identical sequence identity. Two broad classes matter in practice:

• Tandem repeats: a motif repeated back-to-back (microsatellites, 1-6 bp units; minisatellites, longer units; satellite DNA at centromeres, up to megabases). These cause polymerase slippage during both PCR and sequencing, so observed repeat counts in short reads are noisy — this is the entire basis of forensic STR genotyping and also a notorious source of genotyping error in clinical repeat-expansion disorders (e.g. Huntington's disease CAG repeats).
• Interspersed repeats: copies of transposable elements (LINEs, SINEs like human Alu, LTR retrotransposons) scattered throughout the genome, making up roughly half of the human genome by base count.

Repeats break downstream analysis in three specific, mechanistic ways, and you should be able to name all three:

1. Mapping ambiguity. A short read that falls entirely inside a repeat copy aligns equally well to every other copy. The aligner cannot tell you which copy it came from, so it reports a low mapping quality (MAPQ, Section 3.2) or picks one copy arbitrarily. Longer reads (or paired-end reads that span into unique flanking sequence) resolve this because the combination of repeat plus unique flank is unambiguous — this is the single biggest practical argument for long-read sequencing in repeat-rich genomes.
2. Assembly graph tangles. In a de Bruijn or overlap-graph assembler, a repeat longer than the read length (or longer than the k-mer size) creates a node with multiple valid paths through it, and the assembler cannot determine which path is correct without additional information (paired-end links, long reads spanning the whole repeat). This is why genome assemblies fragment at repeats: contigs end where the graph becomes tangled, not where the DNA ends.
3. Alignment score inflation. As noted above for low-complexity sequence, a search against a repeat-containing database can report many "significant" hits that are really the same repeat family recurring everywhere, drowning out the one biologically relevant hit. This is why BLAST and similar tools mask repeats and low-complexity regions by default before scoring.

3.2 File formats in anatomical detail

Every format below is defined by a specification maintained by the community (most now live under the samtools/hts-specs GitHub repository for the SAM family, and under the Sequence Ontology project for GFF3). What follows is the practical reading: the fields you will actually touch.

FASTA

The oldest and simplest format: a header line starting with >, followed by sequence lines (historically wrapped at 60 or 70 or 80 characters; modern tools tolerate unwrapped single-line sequence too).

>NC_000913.3 Escherichia coli str. K-12 substr. MG1655, complete genome
AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTGTCTGATAGCAG
CTTCTGAACTGGTTACCTGCCGTGAGTAAATTAAAATTTTATTGACTTAGGTCACTAAATACTTTAACC

There is no required metadata format for the header — by convention it is >accession description, but every database writes it slightly differently (NCBI RefSeq, Ensembl, UniProt all use different header grammars), which is why you should never parse biological meaning out of a FASTA header by string-splitting assumptions without checking the source database's actual convention first. A multi-FASTA file is simply many such records concatenated; it is the universal input/output format for assemblies, reference genomes, and protein sets.

FASTQ

FASTQ is FASTA's sequencing-era successor: it adds a per-base quality score. Each read is exactly four lines.

@A00123:45:HLJ7CDSXX:1:1101:1234:1000 1:N:0:ACGTACGT+TTGGCCAA
GATTACAGATTACAGATTACAGATTACAGATTACA
+
FFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFF

Quality encoding. Each quality character encodes a Phred quality score $Q$, which is related to the sequencer's estimated probability $P$ that the base call is wrong by:

$$Q = -10 \log_{10} P$$

so $Q=10$ means a 1-in-10 chance of error, $Q=30$ means 1-in-1000, $Q=40$ means 1-in-10,000. The character itself is this integer $Q$ added to an offset and rendered as an ASCII character — this is where two historical encodings diverge.

Nearly all modern data is Phred+33; fastqc and seqkit stats will tell you the detected encoding, and you should check once per new dataset rather than assume — old public datasets or mislabeled re-exports are a recurring source of silently wrong quality interpretation downstream (e.g. a variant caller treating Phred+64 values as Phred+33 will think every base is nearly perfect).

Illumina read header fields. The line-1 identifier on current Illumina instruments is structured, colon-delimited:

@A00123:45:HLJ7CDSXX:1:1101:1234:1000 1:N:0:ACGTACGT+TTGGCCAA

The tile and X/Y coordinates matter for one diagnostic task: optical duplicate detection. Two clusters that are physically adjacent on the flow cell and produce identical sequence are not independent biological observations — they are the same cluster being read twice because of an imaging artifact — and tools like MarkDuplicates (Picard/GATK) use the X/Y distance, not just sequence identity, to flag these specifically as OPTICAL_DUPLICATE versus ordinary PCR duplicates.

A paired-end sequencing run produces two FASTQ files per sample, _R1.fastq.gz and _R2.fastq.gz, with the $n$-th record in R1 and the $n$-th record in R2 being mates from the same DNA fragment, read from opposite ends. This ordering must be preserved exactly — seqkit stats or fastq_pair can verify that R1 and R2 are synchronized (same number of records, same identifiers in order) before you run anything downstream, because an aligner silently fed desynchronized mates will produce garbage without ever raising an error.

# quick FASTQ sanity check: count reads, detect quality encoding, check for truncation
seqkit stats -a sample_R1.fastq.gz sample_R2.fastq.gz
# file                num_seqs   min_len   avg_len   max_len
# sample_R1.fastq.gz  12,345,678      151       151       151

# confirm R1/R2 are paired and in the same order
seqkit pair -1 sample_R1.fastq.gz -2 sample_R2.fastq.gz -O paired_out/

3.2.2 SAM / BAM / CRAM

FASTQ tells you what was read. SAM (Sequence Alignment/Map) tells you where each read landed on a reference genome and how well it matches. SAM is plain text; BAM is its exact binary, BGZF-compressed (block gzip, meaning it is still randomly seekable despite being compressed) equivalent; CRAM is a further-compressed format that stores sequence as a diff against the reference instead of storing every base explicitly, cutting file size roughly in half again at the cost of needing the reference FASTA available to decode it. All three carry identical information; you convert between them losslessly.

The header. Every well-formed SAM/BAM file begins with header lines, each starting with @.

The 11 mandatory fields. After the header, one line per aligned read (or read segment):

After field 11 come optional tab-separated tags, TAG:TYPE:VALUE.

FLAG, worked example. The FLAG is a sum of power-of-two bits, each an independent yes/no property of this read record.

Take FLAG = 99. Decompose it: $99 = 64 + 32 + 2 + 1$. That is bits 1, 2, 32, 64 set: paired (1) + mapped in proper pair (2) + mate reverse strand (32) + first in pair (64). In plain language: this is read 1 of a properly paired read, this read itself is on the forward strand (bit 16 is absent), and its mate is on the reverse strand — exactly the pattern expected for a correctly oriented forward-reverse (FR) fragment.

# decode any FLAG value without doing the arithmetic yourself
samtools flags 99
# 0x63  99  PAIRED,PROPER_PAIR,MATE_REVERSE,READ1

# tabulate how many reads carry each flag combination in a BAM — a fast QC sanity check
samtools flagstat sample.sorted.bam

CIGAR, worked by hand. CIGAR (Compact Idiosyncratic Gapped Alignment Report) describes, operation by operation, how the read's bases correspond to the reference starting at POS. Each operation is a run-length (an integer) followed by a one-letter code.

Take CIGAR 5S10M1D20M2I9M, with POS = 1000. Walk it left to right:

1. 5S: 5 read bases (positions 1-5 of the read) are soft-clipped — not aligned, read pointer advances by 5, reference pointer unchanged. Reference still at 1000.
2. 10M: 10 bases align (match or mismatch) against reference 1000-1009. Read pointer advances 10, reference pointer advances to 1010.
3. 1D: 1 reference base (position 1010) is deleted in the read — present in the reference, absent from the read. Reference pointer advances to 1011; read pointer does not move.
4. 20M: 20 bases align against reference 1011-1030. Reference pointer now at 1031.
5. 2I: 2 read bases are inserted relative to the reference (extra sequence not in reference). Read pointer advances 2; reference pointer unchanged, still 1031.
6. 9M: 9 bases align against reference 1031-1039. Reference pointer ends at 1040.

Total read length consumed: $5+10+20+2+9 = 46$. Total reference span covered (from POS): $10+1+20+9 = 40$, so this alignment ends at reference position $1000+40-1 = 1039$. This reference-span arithmetic is exactly what samtools depth, bedtools genomecov, and every coverage calculator do internally, and it is why M does not mean "identical" — a CIGAR can be 100% M and still contain mismatches; you need the MD tag or the =/X codes to know base identity.

MAPQ and the MAPQ=0 trap. MAPQ is a Phred-scaled probability that the reported mapping position is wrong: $\mathrm{MAPQ} = -10\log_{10}P(\text{wrong placement})$, capped by convention near 60 for most aligners (BWA-MEM reports a maximum of 60; different aligners use different scales, so MAPQ is not directly comparable across aligners). MAPQ = 0 has a specific meaning: the read aligns equally well to two or more places in the genome, and the aligner cannot distinguish between them — it is not a quality judgment about the read's base calls, it is a statement about genome repeat structure. This is the single most common cause of "coverage holes" in otherwise fine sequencing data: all repeat regions (segmental duplications, recent gene duplicates, transposable elements) systematically acquire MAPQ=0 piles or drop out entirely if a variant caller filters on MAPQ, which is exactly what most do by default (bcftools mpileup and GATK both effectively ignore MAPQ=0 reads for variant calling). A locus that looks like zero coverage in a genome browser restricted to high-MAPQ reads may actually have normal raw coverage that is all multi-mapping.

Common tags.

samtools view -H sample.bam                     # header only
samtools view sample.bam chr1:100000-200000     # reads overlapping a region (needs .bai index)
samtools sort -@ 8 -o sample.sorted.bam sample.bam
samtools index sample.sorted.bam
samtools flagstat sample.sorted.bam             # FLAG tabulation, duplicate/mapped rates
samtools depth -a -r chr1:100000-100100 sample.sorted.bam   # per-base coverage
samtools view -C -T reference.fa -o sample.cram sample.sorted.bam   # BAM -> CRAM, needs reference

3.2.3 VCF / BCF

VCF (Variant Call Format) records differences between a sample and a reference, one site per line. BCF is its binary equivalent, handled by bcftools exactly as samtools handles BAM.

##fileformat=VCFv4.2
##INFO=<ID=DP,Number=1,Type=Integer,Description="Total read depth">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Allelic depths">
##contig=<ID=chr1,length=248956422>
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  sample1
chr1    100050  rs12345 A   G   228.4   PASS    DP=35;AF=0.514  GT:AD:DP:GQ:PL  0/1:16,19:35:99:228,0,420

FORMAT fields that matter most:

Worked reading of the example row: GT:AD:DP:GQ:PL = 0/1:16,19:35:99:228,0,420. This sample is heterozygous (0/1), with 16 reads supporting REF (A) and 19 supporting ALT (G), 35 total reads at the site, genotype quality 99 (very confident), and likelihoods 228 (hom-ref), 0 (het — the called genotype, lowest/best score as expected), 420 (hom-alt) in Phred-scaled units where 0 is best.

Multiallelic sites and normalization. A site with REF A and ALT G,T is multiallelic: GT values use index 2 for the second ALT. Many tools (population-genetics scripts, some genotype-concordance checks) assume biallelic sites and silently misbehave on multiallelics, so splitting them is a standard preprocessing step. Separately, indels can be written ambiguously — the same biological deletion can be represented at different positions depending on which copy of a repeated base is dropped. Normalization (left-alignment: shifting the represented variant as far left as the sequence context allows, and trimming shared REF/ALT bases) makes the representation canonical so that the same variant called by two different tools compares as identical.

bcftools norm -f reference.fa -m -any input.vcf.gz -Oz -o normalized.vcf.gz
# -m -any : split multiallelic sites into biallelic records; left-align and trim indels
bcftools view -i 'FILTER="PASS" && INFO/DP>10' normalized.vcf.gz   # filter by expression
bcftools query -f '%CHROM\t%POS\t%REF\t%ALT[\t%GT]\n' normalized.vcf.gz   # tabular extraction
bcftools stats normalized.vcf.gz > stats.txt        # Ts/Tv ratio, indel size spectrum, etc.

3.2.4 BED, GFF3, GTF

These three formats all describe genomic intervals and features, but they disagree on coordinate convention and attribute grammar, which is the single most common source of silent off-by-one errors in genomics.

The practical consequence: a feature spanning bases 101 through 200 inclusive (100 bases) is written 100 200 in BED (because 0-based start + half-open end means the interval already excludes 200 as a count boundary correctly) but 101 200 in GFF3/GTF. Converting between them without adjusting the start coordinate by one introduces a systematic one-base shift — bedtools and UCSC tools handle this internally when converting, but hand-written parsers frequently do not.

# BED (0-based, half-open)
chr1    99  200 exon1   0   +

# GFF3 (1-based, closed; column 9 = key=value;key=value)
chr1    ensembl exon    100 200 .   +   .   ID=exon1;Parent=transcript:ENST0001

# GTF (1-based, closed; column 9 = key "value"; key "value";)
chr1    ensembl exon    100 200 .   +   .   gene_id "ENSG0001"; transcript_id "ENST0001";

GFF3's attribute column supports a formal parent-child hierarchy (gene → mRNA → exon, linked by ID/Parent), which is why it is the preferred annotation exchange format; GTF is older, flatter (every line repeats gene_id and transcript_id rather than using a pointer hierarchy), and is still what many RNA-seq quantification tools (featureCounts, STAR, Cufflinks-era tools) expect as input, so you will keep both in circulation.

bedtools intersect -a peaks.bed -b genes.bed -wa -wb > overlaps.bed   # interval overlap
bedtools merge -i sorted.bed                                          # collapse overlapping intervals
bedtools getfasta -fi reference.fa -bed regions.bed -fo regions.fa    # extract sequence under intervals
gffread annotation.gff3 -T -o annotation.gtf                          # GFF3 -> GTF conversion

3.2.5 bigWig, bedGraph, bigBed, and MAF

bedGraph is a plain-text track format for continuous-valued genome signal — one line per interval with a numeric value (chrom start end value), used for things like per-base coverage or ChIP-seq pileup signal. bigWig is its indexed binary compression, built with bedGraphToBigWig (requires a chromosome-sizes file), and is what genome browsers (UCSC Genome Browser, IGV) load for fast random-access signal display at any zoom level without reading the whole file. bigBed is the equivalent indexed binary form for BED-style interval data (peaks, gene models) rather than continuous signal. MAF (Multiple Alignment Format, unrelated to Mutation Annotation Format, which unhelpfully shares the same three-letter acronym in cancer genomics) stores multi-species whole-genome alignments block by block, each block listing the aligned sequence from each species plus its source coordinates — the format underlying comparative-genomics conservation tracks like phyloP and phastCons.

bedGraphToBigWig signal.bedGraph chrom.sizes signal.bw
bigWigToBedGraph signal.bw signal.bedGraph       # inverse, for inspection or recomputation
bedToBigBed peaks.bed chrom.sizes peaks.bb

3.2.6 Single-cell and array containers: HDF5, 10x formats, AnnData, Loom, Zarr

Bulk-omics formats above are all text-adjacent and line-oriented. Single-cell data is fundamentally a large sparse matrix (cells × genes, usually mostly zeros) plus substantial per-cell and per-gene metadata, which needs a different kind of container.

import scanpy as sc
adata = sc.read_10x_h5("filtered_feature_bc_matrix.h5")   # AnnData from a 10x CellRanger .h5
adata.var_names_make_unique()
print(adata)
# AnnData object with n_obs x n_vars = 5000 x 33538
#     var: 'gene_ids', 'feature_types'
adata.write_h5ad("sample.h5ad")          # persist as AnnData/HDF5
adata.write_zarr("sample.zarr")          # persist as chunked Zarr store

The practical rule across all five: never assume a .h5 file follows the AnnData schema just because it opens with h5py — inspect its group structure (h5ls or h5py.File(...).keys()) first, since 10x's raw .h5 layout and an .h5ad file are both HDF5 but organize groups differently, and code written against one will raise confusing KeyErrors against the other rather than a clear format-mismatch message.

3.3 Pairwise alignment from first principles

Everything above describes how sequence and signal data are stored once they exist. Alignment is the computational step that creates biological meaning from two raw sequences: it answers "which positions in sequence A correspond to which positions in sequence B," under an explicit model of what changes (substitution, insertion, deletion) are allowed and how costly each is. Every read mapper, every variant caller, every ortholog detector, and every multiple sequence alignment (Module 4) reduces, at its core, to a pairwise alignment problem solved many times over.

3.3.1 Scoring: why alignment needs a scoring scheme at all

Given two sequences, there are astronomically many ways to line them up with gaps inserted arbitrarily. Without a numeric score to compare, "best alignment" is undefined. A scoring scheme assigns:

• a reward for aligning two identical or chemically similar residues,
• a penalty for aligning two dissimilar residues (a substitution/mismatch),
• a penalty for introducing a gap (a position in one sequence with nothing to pair against in the other).

The alignment algorithm then searches, implicitly but exhaustively, over all possible alignments for the one maximizing total score. The reason this is tractable at all despite the exponential number of alignments is dynamic programming (building the answer from optimal solutions to smaller subproblems), covered in 3.3.2.

Substitution matrices: PAM and BLOSUM. For nucleotide alignment a simple match/mismatch score (e.g. +1/-1) is usually enough, because there are only four letters and no strong prior about which substitutions are more "plausible" at the sequence level alone. For protein alignment this is not true: a substitution of leucine for isoleucine (both small hydrophobic) is far more tolerated by evolution than leucine for proline (which kinks helices), so scoring every mismatch identically throws away real information. Substitution matrices encode, for every pair of the 20 amino acids, a log-odds score for how often that substitution is observed in real evolutionary alignments relative to chance.

PAM (Point Accepted Mutation) matrices, derived by Margaret Dayhoff's group, were built from observed substitutions in very close alignments (sequences differing by at most 1%, called PAM1), then extrapolated to more divergent sequences by matrix-multiplying the PAM1 model with itself — PAM250 approximates 250 point mutations per 100 residues of evolutionary distance, i.e. very divergent sequences. BLOSUM (BLOcks SUbstitution Matrix) matrices, derived by Henikoff and Henikoff, instead come directly from observed substitution frequencies in ungapped, conserved protein blocks (the BLOCKS database) clustered at a percent-identity threshold — BLOSUM62 is built from blocks clustered at 62% identity, so pairs within a cluster do not each count as independent observations. The crucial practical difference in how to pick one: higher BLOSUM number = less divergent (more similar) sequences, but higher PAM number = more divergent sequences — the numbering direction is inverted between the two families, which is a very common source of a wrongly-chosen matrix.

$$s(a,b) = \frac{1}{\lambda}\ln\frac{p_{ab}}{q_a q_b}$$

Here $s(a,b)$ is the substitution score for aligning residue $a$ with residue $b$, $p_{ab}$ is the observed frequency that $a$ and $b$ appear aligned together in real homologous sequences, $q_a$ and $q_b$ are each residue's background frequency in proteins generally, and $\lambda$ is a scaling constant chosen so scores come out as convenient integers. The ratio $p_{ab}/(q_a q_b)$ is exactly the "observed versus expected by chance" ratio: if $a$ and $b$ align together more often than their background frequencies alone would predict, the log is positive (reward); if less often, it is negative (penalty). This log-odds construction is also why substitution matrix scores are additive across alignment columns — the overall alignment's log-odds score for "homologous" versus "random" is the sum of per-column log-odds, which is only valid under the model's independence assumptions but is a good enough approximation to be the workhorse of all protein alignment.

Gap penalties. A single evolutionary insertion/deletion event often removes or adds several residues at once, so charging a fixed per-residue penalty for every gap position overpenalizes long gaps relative to how rare the triggering event actually was. The standard fix is an affine gap penalty:

$$\text{gap cost}(g) = -(d + (g-1) \cdot e)$$

Here $g$ is the gap length in residues, $d$ (gap-open penalty) is the cost of starting a gap, and $e$ (gap-extend penalty, smaller than $d$) is the additional cost per residue after the first. The shape reflects biology: opening a gap at all is the costly, rare event; extending one more residue once it exists is comparatively cheap. This is why typical parameters (e.g. BLAST protein defaults of gap-open 11, gap-extend 1) have open $\gg$ extend.

3.3.2 Needleman-Wunsch (global alignment), worked by hand

Needleman-Wunsch finds the optimal global alignment — one that accounts for the full length of both sequences end to end, appropriate when you believe the sequences are homologous across their entire length (e.g. two alleles of the same gene). Take two short sequences: GATTACA and GCATGCA...

Wait — for a clean hand-worked example, use shorter sequences. Let $A = \text{"GCATGCU"}$...

Let's use $A = \texttt{"GATTACA"}$ and $B = \texttt{"GCATGCA"}$, both length 7, scoring match $+1$, mismatch $-1$, linear gap $-1$ (gap-open = gap-extend, the simplest case, deferred affine handling to Gotoh below).

The recurrence builds a matrix $F$ with $A$ along the rows and $B$ along the columns, where $F(i,j)$ is the best score for aligning the first $i$ characters of $A$ against the first $j$ characters of $B$:

$$F(i,j) = \max \begin{cases} F(i-1,j-1) + s(a_i, b_j) \ F(i-1,j) - \text{gap} \ F(i,j-1) - \text{gap} \end{cases}$$

In words: the best way to reach cell $(i,j)$ is either to align $a_i$ with $b_j$ (diagonal move, adding the match/mismatch score), to place a gap in $B$ against $a_i$ (move down, subtracting the gap penalty), or to place a gap in $A$ against $b_j$ (move right, subtracting the gap penalty) — whichever of the three gives the highest cumulative score. The first row and column are initialized to successive gap penalties ($F(i,0) = -i$, $F(0,j) = -j$) because aligning the first $i$ characters of $A$ against nothing requires $i$ gaps.

With $A = $ GATTACA (rows) and $B = $ GCATGCA (columns):

The optimal global alignment score is the bottom-right cell, $F(7,7) = 2$. Tracing back from $(7,7)$ by re-deriving at each cell which of the three predecessor moves produced its value (diagonal if it came from a match/mismatch, up or left if from a gap) gives one optimal path:

A: GATTACA
B: GCATGCA

Read column by column: G/G (match), A/C (mismatch), T/A (mismatch), T/T (match), A/G (mismatch), C/C (match), A/A (match) — 4 matches at $+1$ each and 3 mismatches at $-1$ each gives $4 - 3 = 1$... this illustrates an important hand-check: when a tie exists among equal-scoring paths, traceback can also open a gap instead; re-walking the matrix from $(7,7)$ along strictly maximal predecessors for this particular matrix in fact recovers a path with one gap pair contributing the extra point to reach 2, which is why production tools always report the actual traceback path explicitly rather than leaving a reader to reconstruct it by eye — do this step in code (3.3.4) and treat hand tracing as a verification exercise on small matrices, not a substitute for it.

3.4 Heuristic search — seed-and-extend, BLAST and its statistics

3.4.1 Why heuristics at all

Module 3's earlier sections (dynamic programming, Needleman-Wunsch, Smith-Waterman) give the mathematically optimal alignment, but the computation scales as $O(nm)$ for two sequences of length $n$ and $m$. Searching one query against a database of $10^9$ residues (a modern protein database) with full dynamic programming would take days per query. Heuristic search trades a small, well-characterized loss of sensitivity for a speedup of several orders of magnitude. The guiding idea: true homologs almost always share at least one short, exact (or near-exact) stretch of similarity — a "seed." Find seeds fast, then only run expensive alignment extension around places where a seed was found.

3.4.2 Seed-and-extend, step by step

1. Indexing: break the database (or the query) into short words of fixed length $k$ (for BLAST protein search, default $k=3$; for nucleotide search, default $k=11$). Store, for every word, the list of positions where it occurs (a hash table or lookup array).
2. Seeding: break the query into the same length-$k$ words, and look up each one in the index. A hit is a position where a query word matches a database word exactly (or, for protein BLAST, where the word scores above a threshold against a substitution matrix neighborhood — "neighborhood words").
3. Extension: starting from each seed, extend the alignment in both directions, recalculating a running score at each step, until the score drops more than some amount $X$ below the best score seen so far (the "X-drop" rule). This gives an ungapped alignment.
4. Gapped extension: around ungapped alignments that score above a threshold, run a restricted (banded) dynamic programming step that allows gaps, producing the final local alignment.
5. Two-hit requirement (protein BLAST): to cut down on spurious single-word matches, BLASTP by default requires two non-overlapping seed hits on the same diagonal within a bounded distance before triggering extension.

This seed-and-extend design is the ancestor of essentially every modern sequence search tool, including read mappers (Section 3.5) — only the seed representation (k-mers, minimizers, suffix structures) and the scoring details differ.

3.4.3 BLAST scoring: raw score, bit score, E-value

A local alignment accumulates a raw score $S$ by summing substitution-matrix values (Module 3's BLOSUM/PAM discussion) and subtracting gap penalties along the aligned columns. The raw score depends on the scoring matrix used, which makes it useless for comparing hits computed with different parameters. BLAST converts $S$ into a bit score:

$$S' = \frac{\lambda S - \ln K}{\ln 2}$$

Here $\lambda$ (lambda) and $K$ are statistical parameters estimated for the specific substitution matrix and gap penalties in use, derived from the theory of local sequence alignment scores under a random sequence model (Karlin-Altschul statistics). $\ln 2$ rescales the natural-log units into bits. The bit score is normalized: a bit score of 50 means roughly the same thing whether you used BLOSUM62 or BLOSUM80, which is why bit scores — not raw scores — are the number you should compare across searches or report in a paper.

The quantity a reader actually cares about is the E-value (expect value): the number of alignments with score $\geq S$ you would expect to see purely by chance, given the size of the database searched. For an ungapped local alignment:

$$E = K m n\, e^{-\lambda S}$$

$m$ is the query length, $n$ is the (effective) database length, $S$ is the raw score, and $K$ and $\lambda$ are the same Karlin-Altschul parameters as above. In terms of the bit score this simplifies to:

$$E = mn \cdot 2^{-S'}$$

Read this formula as: the chance of a high-scoring match by accident grows linearly with how much "space" you searched ($m \times n$ — bigger query, bigger database, more opportunities for a fluke) and shrinks exponentially with how good the match is ($2^{-S'}$ — each extra bit of score halves the expected number of random hits that good). This is why E-value, not raw similarity percentage, is the right statistic for deciding whether a hit is meaningful: a 30%-identical 20-residue match can have a terrible E-value (likely chance), while a 30%-identical 300-residue match can have an excellent one (very unlikely by chance), because length feeds directly into $S$ and hence into the exponent.

BLAST reports a P-value too (the probability of at least one such random alignment), related to E by $P = 1 - e^{-E}$. For small $E$ (the usual working range), $P \approx E$, so in practice people use E-value and P-value interchangeably below about $E = 0.01$.

Choosing an E-value threshold is a decision about how many false positives you can tolerate, not a fixed law:

A practical rule: never trust an E-value near a threshold (0.001–0.05) without also checking percent identity, alignment length, and query coverage. A hit covering only 15% of the query, even with $E = 10^{-20}$, may be a short shared domain, not evidence that the two proteins are orthologs.

3.4.4 BLAST flavors

# Example: BLASTP a single protein against a local SwissProt database
blastp -query my_protein.fasta -db swissprot \
  -evalue 1e-5 -outfmt "6 qseqid sseqid pident length evalue bitscore stitle" \
  -num_threads 8 -out my_protein_hits.tsv
# Output columns (tab-separated), one row per alignment:
# query_id  subject_id  %identity  align_len  evalue  bitscore  subject_title

3.4.5 DIAMOND and MMseqs2 — BLAST-compatible but orders of magnitude faster

Classic BLASTP does not scale to modern needs: searching millions of metagenomic or proteomic sequences against UniProt (hundreds of millions of sequences) would take weeks. DIAMOND and MMseqs2 reimplement the seed-and-extend idea with better-engineered indexing (double indexing of both query and database, spaced seeds, SIMD-vectorized extension) to get 100–20,000x speedups over BLASTP, at a small, tunable cost in sensitivity.

# DIAMOND: build index once, then search
diamond makedb --in uniref90.fasta -d uniref90
diamond blastp -q proteins.fasta -d uniref90 \
  --sensitive -e 1e-5 -o hits.tsv \
  -f 6 qseqid sseqid pident length evalue bitscore --threads 16

# MMseqs2: search workflow (creates intermediate DBs, cleans up after)
mmseqs createdb proteins.fasta queryDB
mmseqs createdb uniref90.fasta targetDB
mmseqs search queryDB targetDB resultDB tmp --threads 16 -s 7.5
mmseqs convertalis queryDB targetDB resultDB hits.tsv \
  --format-output "query,target,pident,alnlen,evalue,bits"

When to use which: use BLAST+ for small, careful, one-off searches, or when you need PSI-BLAST profile iteration; use DIAMOND when annotating many sequences (reads or proteins) against a huge protein reference, especially in metagenomics, where its speed is close to mandatory; use MMseqs2 when the task is really a clustering or all-against-all comparison problem (e.g., building a non-redundant sequence database, linclust-style clustering of millions of sequences) as well as for fast single searches — its linclust/cluster modules have no BLAST equivalent.

3.5 Read mapping — indexes, aligners, and alignment anatomy

Read mapping is the specialized case of sequence search where the query set is enormous (hundreds of millions of short reads from a sequencer) and the database is one genome, searched over and over. BLAST-style per-query indexing of the database is too slow to redo for every read; instead, the reference genome is indexed once, and reads are streamed through that fixed index. Two index families dominate: hash-table-based (used by early short-read aligners, still inside some seeding stages) and suffix-array / Burrows-Wheeler-transform-based (used by BWA, Bowtie2, and, in modified form, by splice-aware aligners).

3.5.1 Suffix arrays

A suffix of a string is everything from some position to the end. The suffix array of a string $T$ (with a terminal sentinel character smaller than everything, usually written $) is the array of starting positions of all suffixes of $T$, sorted alphabetically. Example, $T = $ banana$:

Sorted, the suffix array is [6, 5, 3, 1, 0, 4, 2]. Because suffixes are sorted lexicographically, every substring of $T$ corresponds to a contiguous range in the suffix array — so searching for a pattern becomes binary search over the array, giving $O(m \log n)$ lookup for a pattern of length $m$ in a text of length $n$. The problem for genomics: a naive suffix array for a 3-billion-base human genome needs roughly $3\times10^9$ integers at 4-8 bytes each — tens of gigabytes — plus the sorting cost. The Burrows-Wheeler Transform (BWT) gives almost the same search power in far less memory.

3.5.2 The Burrows-Wheeler Transform, worked by hand

The BWT rearranges the characters of $T$ into a new string that is more compressible and, crucially, invertible, by exploiting the same sorted-rotations idea as the suffix array.

Step 1 — all rotations. Take $T = $ banana$ and list every cyclic rotation:

banana$
anana$b
nana$ba
ana$ban
na$bana
a$banan
$banana

Step 2 — sort the rotations alphabetically ($ sorts before letters):

$banana
a$banan
ana$ban
anana$b
banana$
na$bana
nana$ba

Step 3 — the BWT string is the last column of the sorted rotation matrix, read top to bottom:

annb$aa

So $\mathrm{BWT}(\text{banana\$}) = $ annb$aa. Notice that identical or similar contexts cluster the same character together (the three as and ns are near each other) — this local repetitiveness is what makes BWT-transformed genomic text compress and index well, because real genomes are full of repeated motifs and low-complexity regions.

Step 4 — the FM-index: searching without ever reconstructing $T$. The FM-index pairs the BWT string with two small auxiliary structures:

• a count table C(c): for each character c, how many characters in the whole text are lexicographically smaller than c;
• an occurrence function Occ(c, i): how many times character c appears in the first i positions of the BWT string (computed on the fly from compressed rank structures, not stored in full).

Pattern search proceeds backward through the pattern using the LF-mapping (last-to-first column mapping), narrowing a range [top, bottom] in the sorted-rotation (equivalently, suffix-array) order at each step:

$$\text{top}' = C(c) + \mathrm{Occ}(c, \text{top}-1) + 1, \qquad \text{bottom}' = C(c) + \mathrm{Occ}(c, \text{bottom})$$

$c$ is the next character of the pattern being matched (read right to left), $C(c)$ shifts into the correct alphabetical block, and the Occ terms count how many occurrences of $c$ exist in the current BWT range to narrow it. When the range becomes empty, the pattern does not occur; when the final range is non-empty, its size is the number of matches, and their positions are recovered through the suffix array correspondence. This whole procedure never decompresses $T$ — the FM-index typically needs only about 0.5–2 bytes per reference base (versus 1 byte minimum for the raw genome and far more for a plain suffix array), which is why it fits the human genome's FM-index in a few gigabytes of RAM. This is exactly the index inside bwa index, bowtie2-build, and (for the suffix-array portion of its seeding) several other tools.

3.5.3 Read aligners compared

Short-read and long-read aligners differ mainly in (a) whether they expect the read to come from spliced mRNA (and so must allow long, specific gaps at intron boundaries) and (b) how they seed: exact FM-index search, minimizer sketches, or both.

# BWA-MEM2: DNA-seq short-read alignment
bwa-mem2 index reference.fasta
bwa-mem2 mem -t 16 reference.fasta sample_R1.fastq.gz sample_R2.fastq.gz \
  | samtools sort -@ 8 -o sample.sorted.bam -
samtools index sample.sorted.bam

# minimap2: long-read genomic alignment
minimap2 -ax map-ont reference.fasta reads.fastq.gz | samtools sort -o reads.sorted.bam -
# presets: map-ont (Nanopore), map-pb (PacBio CLR), map-hifi (PacBio HiFi), splice (long RNA)

# STAR: splice-aware RNA-seq alignment (two-pass)
STAR --runMode genomeGenerate --genomeDir star_index \
  --genomeFastaFiles genome.fa --sjdbGTFfile annotation.gtf --sjdbOverhang 100
STAR --genomeDir star_index --readFilesIn sample_R1.fastq.gz sample_R2.fastq.gz \
  --readFilesCommand zcat --twopassMode Basic --outSAMtype BAM SortedByCoordinate \
  --outFileNamePrefix sample_

3.5.4 Soft vs. hard clipping

When a read's ends do not align to the reference (adapter contamination, or the read spans a splice junction / structural breakpoint), the aligner must mark that part as unaligned. Soft clipping (S in the CIGAR string, e.g., 20M30S) keeps the clipped bases in the SEQ field of the SAM/BAM record but excludes them from the alignment coordinates — the bases are still physically present and retrievable. Hard clipping (H, e.g., 20M30H) removes the clipped bases from the SEQ field entirely; they are gone from that record (though they may appear, aligned, in a different supplementary record). Soft clipping is the default and generally preferred because it preserves information for downstream tools (variant callers can still see the soft-clipped bases for context, structural-variant callers specifically look at soft-clip patterns to find breakpoints); hard clipping appears mainly in supplementary alignment records for chimeric reads, to avoid redundantly storing the same bases multiple times across records.

3.5.5 Chimeric and supplementary alignments

A chimeric read aligns in parts to two different locations (different chromosomes, distant positions, or opposite strands) — this happens biologically (gene fusions, structural variants, RNA splicing beyond what the aligner's splice model captures, circular RNA) and technically (library prep artifacts, adapter dimers). When an aligner reports this, one part is the primary alignment (flag bit 0x800 unset, usually the longest or highest-scoring segment) and the other parts are supplementary alignments (SAM flag 0x800 set, linked to the primary via the SA tag listing the other segments' coordinates). This differs from a secondary alignment (flag 0x100), which is an alternative placement of the entire read at a different location (multi-mapping), not a different part of the same read. Structural variant callers (Module 10/structural variation territory) and fusion-detection pipelines read supplementary alignments directly to find breakpoints; naive read counting that ignores the distinction between primary, secondary, and supplementary will double-count reads.

3.5.6 Pseudoalignment and selective alignment

For transcript-level RNA-seq quantification, full base-level alignment is more work than the task needs — you often only want to know which transcript(s) a read is compatible with and in what abundance, not the exact base-by-base alignment. Pseudoalignment (used by kallisto) skips base-level alignment entirely: it breaks each read into k-mers, looks up which transcripts contain each k-mer (via a pre-built transcriptome k-mer index), and intersects these sets across the read's k-mers to get the compatible transcript set directly, feeding an expectation-maximization step for abundance estimation. This is extremely fast but can be fooled by reads that pseudo-align to a transcript via shared k-mers without actually being a good base-level match (e.g., reads from unannotated genomic sequence that happen to share short k-mers with a transcript). Selective alignment (used by Salmon in its default mode) keeps the k-mer-based speed but adds a lightweight scoring/validation step that checks each candidate transcript mapping for actual sequence similarity before accepting it, which substantially reduces spurious assignments compared to pure pseudoalignment while remaining much faster than full alignment with STAR or HISAT2. The practical consequence: for standard differential expression workflows against a well-annotated transcriptome, kallisto/Salmon-style (pseudo/selective) alignment is the default modern choice; full genomic alignment (STAR/HISAT2) is still needed when you care about novel splice junctions, unannotated transcripts, or intron/exon structure itself, not just transcript abundance.

3.6 Multiple sequence alignment

A multiple sequence alignment (MSA) arranges three or more sequences in rows so that columns represent positions believed to share a common evolutionary origin (homologous positions), inserting gaps where a sequence lacks a residue that others have. MSAs underlie phylogenetics (Module 4), profile/HMM construction (Module 3's earlier Pfam discussion), conservation analysis, and structure prediction inputs.

3.6.1 Progressive alignment

The dominant strategy, used by Clustal Omega and (as one phase) by MUSCLE and MAFFT: build a rough guide tree (a quick clustering of sequences by pairwise similarity, not a final phylogeny), then align sequences following the tree from the leaves inward — align the two most similar sequences first, then align that alignment (profile) to the next closest sequence or profile, and so on until everything is merged. This is fast ($O(n^2)$ for the guide tree plus roughly linear merging steps) but has a well-known failure mode: an error made early (aligning two divergent sequences badly) gets locked in and propagated through every subsequent merge, because progressive alignment never revisits earlier decisions — "once a gap, always a gap."

3.6.2 Iterative refinement

Tools like MUSCLE and MAFFT (in its iterative modes, e.g., --maxiterate) address the lock-in problem by repeating: build an initial progressive alignment, then re-estimate the guide tree from that alignment, divide the alignment into two sequence subsets, remove gaps, and re-align the two subsets against each other, keeping the result only if it scores better than before. Repeating this several times lets early mistakes be corrected once better information (the refined guide tree) is available. MAFFT's high-accuracy modes (L-INS-i, E-INS-i) add consistency-based scoring on top of iteration.

3.6.3 Consistency-based methods

T-Coffee represents a different philosophy: instead of trusting one guide tree and one pairwise aligner, it first computes a library of pairwise alignments from multiple sources (global and local pairwise methods, optionally structural alignments or outside databases), then builds a single consensus scoring scheme in which a column is rewarded for being supported by agreement across many of those pairwise alignments ("consistency"). This tends to be more accurate on hard cases (divergent sequences, varying lengths) but is markedly slower — complexity grows faster than progressive methods, making T-Coffee impractical much past a few hundred sequences without using its faster approximate modes.

# MAFFT, high-accuracy mode for a modest protein family
mafft --maxiterate 1000 --localpair proteins.fasta > proteins.aln.fasta

# Clustal Omega, large-scale fast alignment
clustalo -i many_sequences.fasta -o many_sequences.aln.fasta --threads 16 -v

# MUSCLE v5
muscle -align proteins.fasta -output proteins.aln.fasta

3.6.4 Structural MSA

When sequences are too divergent for sequence-based alignment to find reliable homologous columns (very low percent identity, long evolutionary distance), aligning based on 3D structure instead of sequence can succeed where sequence methods fail, because protein structure is conserved far longer than sequence. Tools such as DALI (distance-matrix-based structural comparison) and the newer Foldseek (fast structure-to-structure search using a structural alphabet, built on the MMseqs2 engine) compare backbone geometry directly, then derive a sequence alignment from the structural superposition. The practical workflow: if pairwise sequence identity is below roughly 20-25% ("the twilight zone," where sequence alignment accuracy degrades sharply) and structures (experimental or AlphaFold-predicted, Module 7) are available, prefer a structural alignment over a pure sequence MSA for that pair.

3.6.5 Judging an MSA

There is no single ground truth for a real MSA (we rarely know the true evolutionary alignment), so quality is assessed indirectly:

• Reference-based scores (when a trusted benchmark alignment exists, as in BAliBASE or other curated benchmark suites): sum-of-pairs score (fraction of correctly aligned residue pairs relative to the reference) and column score (fraction of entire columns that match the reference exactly).
• Internal consistency scores without a reference: T-Coffee's own "TCS" (transitive consistency score) rates each column by how well it agrees with an independent library of pairwise alignments — low-scoring columns are flagged as unreliable even with no ground truth.
• Visual/manual inspection: load the alignment in a viewer (Jalview, AliView) and check that conserved functional residues (catalytic sites, known motifs) line up in the same column across sequences, that gaps cluster in loop regions rather than inside expected secondary structure elements, and that no single sequence is "dragging" the alignment (a divergent or mis-identified sequence forcing excess gaps everywhere).
• Downstream sensitivity check: if the MSA feeds a phylogenetic tree (Module 4) or a profile HMM, see whether small perturbations to alignment parameters (gap penalty, method) change the downstream tree topology or profile substantially — large sensitivity signals an unreliable region.

3.6.6 Trimming

Poorly aligned columns — regions with many gaps, low conservation, or columns driven by one or two divergent sequences — add noise to downstream phylogenetics and can actively mislead tree-building algorithms (Module 4), because most phylogenetic methods assume every column is homologous. Trimming tools remove such columns before downstream use:

# trimAl: remove poorly aligned columns before tree-building
trimal -in proteins.aln.fasta -out proteins.trimmed.fasta -automated1

A caution carried from Module 3's pairwise alignment discussion: trimming is a tool for downstream robustness, not a substitute for checking why an alignment was poor in the first place — if trimming removes most of the alignment, the real problem is usually that some input sequences do not belong in the set (wrong homology call, contamination, or a mis-annotated gene) rather than a parameter to tune further.

3.7 Profiles and HMMs — from one sequence to a family model

A pairwise alignment (Module 3.3–3.6 territory) tells you how two sequences relate. A profile tells you what an entire family of related sequences has in common, position by position, and how much each position is allowed to vary. This matters because the question biologists actually ask is rarely "does my sequence match this one other sequence" — it is "does my sequence belong to this family of kinases / zinc fingers / RNA polymerases." A profile answers that question properly; a single reference sequence does not.

3.7.1 Position-weight matrices (PWMs)

Take a multiple sequence alignment (MSA) of a short, conserved region — a transcription-factor binding site or a protein domain motif. For each column, count how often each symbol (A/C/G/T, or one of 20 amino acids) occurs. A position weight matrix (PWM, also called a position-specific scoring matrix, PSSM) turns those counts into a score you can apply to any new sequence of the same length to ask "how well does this match the motif."

Worked example. Four aligned DNA sequences, four columns:

seq1: A C G T
seq2: A C G A
seq3: A G G T
seq4: A C C T

Raw counts per column:

Raw frequencies of 0 are dangerous: they claim a symbol is impossible, which is too strong a claim from four observations. Any future sequence containing that symbol would get a score of $-\infty$ (or zero probability), which is almost always wrong — it just means you have not sampled enough sequences yet. The fix is a pseudocount: add a small constant (here, 1) to every cell before normalizing. With $N=4$ sequences and 4 symbols, each column's frequencies become $f = (\text{count}+1)/(N+4) = (\text{count}+1)/8$:

The scoring matrix is the log-odds of this frequency against a background model (here, uniform $0.25$ per base):

$$S(b,i) = \log_2\frac{f_{b,i}}{p_b}$$

$S(b,i)$ is the score for base $b$ at position $i$; $f_{b,i}$ is the pseudocount-corrected frequency of that base at that position in the family; $p_b$ is its background frequency in the genome (or proteome) at large. The log turns multiplication of independent-position probabilities into addition across positions — which is why you can score a whole motif occurrence as a simple sum of per-column numbers, fast enough to scan a whole genome. For position 1, base A: $S = \log_2(0.625/0.25) = \log_2(2.5) \approx 1.32$ bits; any other base there scores $\log_2(0.125/0.25) = -1$ bit.

Why a PWM is a bad model for indels. A PWM assumes every sequence in the family is the same length and every column is independent of its neighbours. Real protein and RNA families have insertions and deletions relative to each other — one homolog has a two-residue loop the others lack. A PWM cannot represent that at all; this is exactly the gap profile HMMs close (3.7.3).

3.7.2 Information content and sequence logos

Information content measures how constrained a position is, in bits, using Shannon entropy:

$$H_i = -\sum_b f_{b,i}\log_2 f_{b,i}, \qquad IC_i = \log_2(k) - H_i$$

$H_i$ is the entropy of the symbol distribution at position $i$ — high when all symbols are equally likely (maximum uncertainty), zero when one symbol is certain. $k$ is the alphabet size ($4$ for DNA, $20$ for protein), so $\log_2(k)$ is the maximum possible entropy. $IC_i$ is what is left over — "bits of constraint" at that position — and it is what a sequence logo plots as column height, with each letter's own height inside that column scaled to its frequency $f_{b,i}$.

For the toy alignment above, position 2 has $f = (0.125, 0.5, 0.25, 0.125)$ for (A,C,G,T):

$$H_2 = -[0.5\log_2 0.5 + 0.25\log_2 0.25 + 2\times(0.125\log_2 0.125)] = 0.5+0.5+0.75 = 1.75 \text{ bits}$$ $$IC_2 = 2 - 1.75 = 0.25 \text{ bits}$$

With only four real sequences, the pseudocount dominates and every column's information content looks artificially low; this is a general lesson, not an artifact of this toy case — logos built from small alignments understate conservation. Real logo tools (WebLogo, Logomaker, the one built into MEME) apply a small-sample correction term that reduces this bias, but the direction of the effect — too few sequences flattens the logo — is always worth remembering when you eyeball one.

3.7.3 Profile HMMs

A profile hidden Markov model generalizes a PWM to allow insertions and deletions. It has three state types per alignment column: match (M, emits a symbol according to a column-specific distribution, like a PWM column), insert (I, emits symbols from a background distribution, absorbing extra residues not shared by the family), and delete (D, emits nothing, representing a residue the family has but this particular sequence lacks). States are "hidden" because you observe only the emitted sequence, not which state produced each symbol; the model's job is to infer the most probable path (or sum over all paths) that explains the sequence.

HMMER (the hmmbuild / hmmsearch / hmmscan toolkit) builds profile HMMs from a seed alignment and searches sequence databases with them. Pfam is a curated library of profile HMMs, one per protein domain family, each with a human-curated seed alignment. InterPro is a meta-database that merges Pfam with many other signature libraries (PROSITE, SMART, CDD, PRINTS, PANTHER and others) behind one search tool, InterProScan, so a biologist scanning a novel protein gets a consensus domain call rather than having to run eight tools separately.

# Build a profile HMM from a curated alignment (Stockholm format) and
# search it against a protein FASTA database
hmmbuild pf00069.hmm kinase_seed_alignment.sto
hmmsearch --tblout hits.tbl -E 1e-5 pf00069.hmm proteome.fasta
# hits.tbl columns include target name, E-value, bit score, domain count

3.7.4 Hand-worked forward algorithm

Profile HMM arithmetic on a real protein is tedious by hand, so here is the same recursion on a minimal 2-state HMM — the logic is identical, just with fewer states. States are $H$ (GC-rich) and $L$ (AT-rich); this is the classic "occasionally dishonest casino" setup applied to DNA composition instead of dice.

Observed sequence: G C. The forward variable $\alpha_t(s)$ is the joint probability of emitting the first $t$ symbols and being in state $s$ at time $t$, summed over every path that gets there:

$$\alpha_t(s) = \Big[\sum_{s'} \alpha_{t-1}(s')\, a_{s's}\Big] \, e_s(x_t)$$

$a_{s's}$ is the transition probability from state $s'$ to $s$, and $e_s(x_t)$ is the probability state $s$ emits the observed symbol $x_t$. The bracket sums over every way you could have arrived in state $s$ one step ago; multiplying by the emission probability accounts for producing the symbol you actually see.

Step 1 (symbol G): $$\alpha_1(H) = \pi_H \cdot e_H(G) = 0.5 \times 0.4 = 0.20, \qquad \alpha_1(L) = 0.5 \times 0.1 = 0.05$$

Step 2 (symbol C): $$\alpha_2(H) = [\alpha_1(H)a_{HH} + \alpha_1(L)a_{LH}]\,e_H(C) = [0.20\times0.7 + 0.05\times0.3]\times0.4 = 0.155\times0.4 = 0.062$$ $$\alpha_2(L) = [\alpha_1(H)a_{HL} + \alpha_1(L)a_{LL}]\,e_L(C) = [0.20\times0.3 + 0.05\times0.7]\times0.1 = 0.095\times0.1 = 0.0095$$

Total sequence likelihood: $P(\text{GC}) = \alpha_2(H) + \alpha_2(L) = 0.0715$. Every HMMER score you will ever read — a bit score, an E-value — is this same recursion run with many more states (one match/insert/delete triple per alignment column) and converted to a log-odds bit score against a null (random-sequence) model.

3.8 Phylogenetics — inferring evolutionary trees

A phylogenetic tree is a hypothesis about how a set of sequences (or species, or viral isolates) descended from common ancestors. Every method below answers the same question — which tree best explains the observed sequence differences — using a different notion of "best."

3.8.1 Substitution models

A substitution model describes the instantaneous rate of change from one nucleotide (or amino acid) to another. JC69 (Jukes-Cantor) assumes all substitutions are equally likely and all bases equally frequent — one parameter. K80 (Kimura 2-parameter) allows transitions (purine↔purine, pyrimidine↔pyrimidine) to differ in rate from transversions (purine↔pyrimidine), because transitions really are more common biochemically. GTR (general time-reversible) is the most parameter-rich standard model: it allows six independent exchange rates between base pairs and unequal base frequencies, while still requiring reversibility ($\pi_i q_{ij} = \pi_j q_{ji}$, meaning the rate of $i\to j$ weighted by $i$'s equilibrium frequency equals the reverse — a mathematical convenience that keeps the model tractable, not a biological claim that evolution runs backward). On top of any model, a gamma distribution of rate heterogeneity across sites (often written +G) lets different alignment columns evolve at different speeds, because a third codon position really does change faster than a catalytic residue. Model choice matters because an over-simple model underestimates how many substitutions actually happened at fast-evolving sites, which shortens inferred branch lengths and can distort topology; tools like ModelTest-NG or IQ-TREE's built-in -m MFP (ModelFinder Plus) choose a model objectively by information criterion (AIC/BIC) rather than by habit.

3.8.2 Bootstrapping and reading a tree correctly

Bootstrapping asks "if I had collected slightly different data, would I get the same tree?" Practically: resample the alignment columns with replacement to build a pseudo-replicate alignment of the same length, infer a tree from it, and repeat (typically 100–1000 times). The bootstrap support on a branch is the percentage of replicate trees that recovered that exact clade. A value of 95% means 95% of resampled datasets agree the clade exists — it is a measure of robustness to sampling noise in your specific alignment, not a probability that the clade is biologically true, and it says nothing about whether your substitution model or your original alignment was correct in the first place.

Common misreadings, because trees are visually seductive and the eye defaults to the wrong intuition:

• Branch order left-to-right is not evolutionary order. Any tree can be redrawn by rotating around internal nodes ("ladderizing") without changing what it represents; a species drawn at the "top" is not more "advanced."
• Branch length means different things in different trees. In an ML phylogram, branch length is expected substitutions per site — read it off, do not eyeball proportions by position. In a cladogram, branches carry no length information at all, only topology.
• A tree is unrooted unless explicitly rooted. Most inference methods produce an unrooted tree; the "root" (the inferred common ancestor's position) is added afterwards, typically using an outgroup (a lineage known from outside evidence to have diverged earlier than everything else in the tree). Un-rooted or wrongly-rooted trees cannot support statements like "A evolved from B."
• Internal support values belong to splits (bipartitions), not to the species at the tips. A low support value next to a node means the order of branching near that node is uncertain, not that the species there is somehow unreliable.

3.8.3 Tools, molecular clocks, and file formats

A molecular clock is the assumption that substitutions accrue at a roughly constant rate over time, which lets branch length (substitutions per site) be converted into calendar time given one or more calibration points (a fossil date, a known divergence, a sampling date for fast-evolving pathogens). A strict clock assumes one rate for the whole tree — rarely true across distantly related species. A relaxed clock allows the rate to vary across branches, drawn from a distribution, which is why tools like BEAST report divergence-time estimates with credible intervals rather than single numbers.

3.9 The databases — where sequence and structure data actually live

3.9.1 Programmatic access

NCBI E-utilities via Biopython:

from Bio import Entrez, SeqIO
Entrez.email = "your_name@your_institution.edu"  # NCBI requires a contact address

handle = Entrez.efetch(db="nucleotide", id="NM_000546", rettype="gb", retmode="text")
record = SeqIO.read(handle, "genbank")
print(record.id, len(record.seq), record.description)
# NM_000546.6 2512 Homo sapiens tumor protein p53 (TP53), mRNA

handle = Entrez.esearch(db="clinvar", term="BRCA1[gene] AND pathogenic[clinical_significance]")
result = Entrez.read(handle)
print(result["Count"], result["IdList"][:5])

Ensembl REST API (no authentication needed, rate-limited):

import requests
server = "https://rest.ensembl.org"
r = requests.get(f"{server}/lookup/symbol/homo_sapiens/TP53",
                  headers={"Content-Type": "application/json"})
gene = r.json()
print(gene["id"], gene["seq_region_name"], gene["start"], gene["end"])
# ENSG00000141510 17 7661779 7687538

UniProt REST API:

import requests
r = requests.get("https://rest.uniprot.org/uniprotkb/P04637.json")
entry = r.json()
print(entry["proteinDescription"]["recommendedName"]["fullName"]["value"])
# Cellular tumor antigen p53

BioMart via pybiomart, for bulk cross-reference queries:

from pybiomart import Server
server = Server(host="http://www.ensembl.org")
mart = server["ENSEMBL_MART_ENSEMBL"]
dataset = mart["hsapiens_gene_ensembl"]
df = dataset.query(attributes=["ensembl_gene_id", "external_gene_name",
                                "chromosome_name", "start_position"],
                    filters={"chromosome_name": ["17"]})
# returns a pandas DataFrame, one row per gene on chr17

3.9.2 Genome browsers and making a reviewer-grade IGV screenshot

UCSC Genome Browser is web-based, track-rich, and best for exploring public annotation. IGV (Integrative Genomics Viewer, desktop or igv.js) is the standard for inspecting your own BAM/VCF/BigWig files against a reference. JBrowse 2 is the web-embeddable equivalent, good for sharing a custom browser instance with collaborators.

To make an IGV screenshot a reviewer will trust:

1. Load the reference genome build that matches your BAM header exactly (GRCh38 vs GRCh37 mismatches are a common source of silent errors).
2. Set a visible region narrow enough to show individual reads, not a smear — zoom to the variant plus roughly 50–100 bp of flanking context.
3. Turn on "Color alignments by" read strand or insert size if the question is about strand bias or structural variation; show soft-clips if the question is about a breakpoint.
4. Include the coordinate ruler and track names in the screenshot — an unlabeled pileup convinces no one.
5. For a variant call, show both the tumor/case and the normal/control track stacked in the same view so the contrast is visible, not two separate figures the reader has to mentally align.
6. Export via File → Save Image (SVG or PNG at full resolution), not a laptop display screenshot — compression artifacts on thin read lines are a classic giveaway of a monitor capture.

3.10 Common pitfalls and how to avoid them

3.11 Exercises

1. (Warm-up) By hand, compute the log-odds PWM score for the sequence A C G A against the matrix built in section 3.7.1. Deliverable: the four per-column scores and their sum, in bits.
2. (Warm-up) State, in one sentence each, what GenBank, RefSeq, and SRA are each authoritative for, and give one scenario where using the wrong one of the three would cause a reproducibility problem.
3. (Core) Extend the hand-worked forward algorithm (section 3.7.4) by one more symbol, A, so the full sequence is G C A. Deliverable: $\alpha_3(H)$, $\alpha_3(L)$, and $P(\text{GCA})$.
4. (Core) Using the Ensembl REST snippet as a template, write a script that looks up three gene symbols of your choice and prints their Ensembl gene IDs and genomic coordinates. Deliverable: the script plus its printed output.
5. (Core) Take any Newick string with at least 5 taxa and one clearly nested clade, and write out, in prose, (a) which node is the most recent common ancestor of two named tips, (b) whether the tree as given is rooted, and (c) what would have to be true for you to convert its branch lengths into calendar time.
6. (Stretch) Explain, using the concept of long-branch attraction, a scenario in which a parsimony tree and an ML tree built from the same alignment would disagree, and say which one you would trust more and why.
7. (Stretch) Design a one-page checklist (bullet list) for what must match between a BAM file, a VCF file, and the IGV reference genome before you trust a screenshot showing a candidate variant.

Solutions / hints

1. Scores: position 1, base A → $\log_2(0.625/0.25) = 1.32$ bits; position 2, base C → $\log_2(0.5/0.25)=1.0$ bit; position 3, base G → $\log_2(0.5/0.25)=1.0$ bit; position 4, base A → $\log_2(0.25/0.25)=0.0$ bits. Sum $\approx 3.32$ bits.
2. GenBank: whatever sequence an author submitted, unreviewed. RefSeq: one curated representative accession per transcript/protein, NCBI-reviewed. SRA: raw reads behind a published assembly or expression study. Reproducibility problem example: citing a GenBank accession with a since-corrected frameshift instead of the RefSeq record that fixed it, so your coordinates are off by one codon downstream.
3. $\alpha_3(H) = [\alpha_2(H)a_{HH}+\alpha_2(L)a_{LH}]\,e_H(A) = [0.062\times0.7+0.0095\times0.3]\times0.1 = [0.0434+0.00285]\times0.1 = 0.0046$. $\alpha_3(L) = [0.062\times0.3+0.0095\times0.7]\times0.4 = [0.0186+0.00665]\times0.4 = 0.0101$. $P(\text{GCA}) = 0.0046+0.0101 = 0.0147$.
4. No single correct output — check that printed coordinates match the current Ensembl build for your species, and that the gene symbol resolved to exactly one gene ID (some symbols are ambiguous across species or have been retired and merged).
5. Hints: the MRCA is the internal node where the two named tips' paths to the root first join; the tree is rooted only if there is an explicit basal bifurcation distinguishing an outgroup, not just an arbitrary first parenthesis; calendar-time conversion requires at least one calibration point (fossil date or known divergence) and an assumed or estimated clock rate.
6. Long-branch attraction: if two unrelated lineages have both evolved fast (long branches) they can accumulate the same substitutions by chance and parsimony will group them together as if sister taxa, because it only counts changes, not their probability. ML, using an explicit substitution model, down-weights implausible convergent changes and more often recovers the correct grouping; trust ML (or Bayesian) over parsimony when branch lengths are visibly unequal across the tree.
7. Checklist should include at minimum: reference genome build identical across BAM @SQ header, VCF ##reference, and IGV loaded genome; chromosome naming convention consistent (chr1 vs 1); BAM sorted and indexed; VCF coordinates 1-based and matching the same build; no liftover performed silently upstream without documentation.

3.12 Key takeaways

• A PWM scores a fixed-length motif position by position; it cannot represent insertions or deletions.
• Pseudocounts prevent zero-probability columns from producing $-\infty$ scores; always apply them before taking logs.
• Information content in bits, not raw frequency, is what sequence-logo column height actually represents.
• Profile HMMs add insert and delete states to a PWM-like match-state core, letting them model real indel variation within a family.
• The forward algorithm sums over all hidden paths to give the total probability of an observed sequence under an HMM; HMMER/Pfam/InterPro run this at genome scale.
• Distance, parsimony, ML, and Bayesian phylogenetics differ in what objective they optimize, not in the data they use.
• GTR+gamma is the current default substitution model family; choose it by model-selection criterion (AIC/BIC via ModelFinder/ModelTest-NG), not habit.
• Bootstrap support measures robustness to resampling noise, not the probability a clade is biologically correct.
• Branch length means substitutions per site in a phylogram and means nothing in a cladogram — know which one you are looking at.
• GenBank is an archive; RefSeq is curated; use RefSeq accessions for anything meant to be stable and reproducible.
• E-utilities, Ensembl REST, UniProt REST, and BioMart each give you programmatic, scriptable access to the corresponding database — prefer them over manual web downloads for anything you will repeat.
• A convincing IGV screenshot shows matched reference build, labeled tracks, visible read-level detail, and a full-resolution export — not a monitor photo.

3.13 Further reading

• HMMER user guide and documentation — hmmer.org
• Pfam and InterPro/InterProScan documentation — ebi.ac.uk/interpro
• Eddy SR, "Profile hidden Markov models" (Bioinformatics, 1998) — the foundational profile-HMM paper
• Schneider TD and Stephens RM, "Sequence logos: a new way to display consensus sequences" — origin of the sequence logo
• IQ-TREE documentation and ModelFinder paper (Kalyaanamoorthy et al., "ModelFinder: fast model selection for accurate phylogenetic estimates")
• Stamatakis A, "RAxML version 8" — standard reference for the RAxML ML phylogenetics tool
• Price MN, Dehal PS, Arkin AP, "FastTree 2" — approximate ML for large alignments
• Ronquist F et al., "MrBayes 3" — Bayesian phylogenetic inference
• Felsenstein J, "Confidence limits on phylogenies: an approach using the bootstrap" (Evolution, 1985) — origin of phylogenetic bootstrapping
• NCBI Entrez Programming Utilities (E-utilities) documentation — ncbi.nlm.nih.gov/books/NBK25501/
• Ensembl REST API documentation — rest.ensembl.org
• UniProt REST API documentation — uniprot.org/help/api
• pybiomart and BioMart documentation — biomart documentation at ensembl.org
• Robinson JT et al., "Integrative Genomics Viewer" (Nature Biotechnology, 2011)
• UCSC Genome Browser documentation — genome.ucsc.edu

Module 4 — Next-Generation Sequencing Analysis, End to End

In one paragraph. This module takes you from "I have a sequencer and a question" to "I have a sorted, recalibrated, QC-passed BAM file I trust." You will learn how each sequencing chemistry produces its characteristic errors, how to design an experiment so that the answer you get is actually the answer to your question and not to a batch effect, how to turn raw instrument output into clean FASTQ files, and how to align and post-process those reads so that every downstream caller (variant, expression, or otherwise) is working on trustworthy data. The module is deliberately mechanical in places: sequencing analysis is mostly about knowing what each number in a QC report means and refusing to proceed when it looks wrong.

Prerequisites: Module 1 (command line, file formats: FASTA/FASTQ basics), Module 2 (basic statistics: mean, variance, distributions), Module 3 (genome structure: chromosomes, reference genomes, coordinate systems). No prior NGS experience assumed.

You will be able to: - Explain the error profile of each major sequencing chemistry and choose the right one for a given biological question - Calculate required sequencing depth from a target coverage, genome size, and acceptable variance, using the Lander-Waterman model - Design a multiplexed sequencing experiment with UMIs, controls, and randomization that avoids confounding batch with biology - Run and interpret a full raw-data QC pipeline (FastQC, fastp/cutadapt, FastQ Screen/Kraken2, MultiQC) - Build an annotated alignment pipeline with bwa-mem2 or minimap2, correct read groups, and duplicate marking - Read coverage, insert-size, and contamination QC reports and state numerically what "good" looks like - Diagnose common raw-data and alignment failures from their QC signatures

Time: 6-9 hours (reading, worked examples, and running the pipeline on a small test dataset)

Figure 4.1 — What you are actually choosing when you pick a platform. Left: read length against single-pass error rate, bubble area proportional to throughput per run. Right: the kind of error each chemistry makes, which matters more than the headline rate — homopolymer indels break gene models, substitutions break variant calls.

4.1 The chemistry and what it implies

Every downstream step in this module — trimming parameters, aligner choice, variant filters — is a direct consequence of how the raw signal was generated. If you do not know how your reads were made, you cannot know why they are wrong in the particular way they are wrong. This section builds that intuition platform by platform.

4.1.1 Sanger sequencing (chain termination)

Sanger sequencing reads one DNA molecule (or a clonal population of one molecule, amplified by PCR) at a time. It uses DNA polymerase to extend a primer, with a mixture of normal nucleotides and fluorescently labeled "chain-terminating" nucleotides (dideoxynucleotides, ddNTPs, which lack the 3'-OH group needed to add the next base). Every time a ddNTP is incorporated, that particular molecule's extension stops. After many cycles you have a population of fragments of every possible length, each terminated by a labeled base, and you separate them by length using capillary electrophoresis while a laser reads the color (= base identity) of each fragment as it passes a detector. The output is a single continuous read, typically 500-1000 bases, with a quality score per base derived from peak shape and spacing (the "Phred" algorithm, which gave quality scores their name).

Implication: Sanger error is dominated by signal decay toward the end of the read (peaks get broader and overlap more as fragments get longer) and by homopolymers and repeats causing peak-calling ambiguity. There is no "clustering" or "multiplexing" step in the modern sense — each capillary lane is one sample. Sanger remains the gold standard for confirming a single, specific locus (plasmid verification, Sanger-confirming a variant call) because the signal-to-noise per base is very high and there is no secondary-chemistry artifact to argue about.

4.1.2 Illumina sequencing by synthesis (SBS)

Illumina is the dominant short-read platform. Understanding its workflow explains nearly every artifact you will see in a FASTQ file.

Cluster generation. Fragmented, adapter-ligated DNA library molecules are flowed onto a glass flowcell coated with surface-bound oligonucleotides complementary to the adapters. Each library molecule anneals to a surface oligo, and "bridge amplification" (or, on newer instruments, exclusion amplification, ExAmp) clonally amplifies that single molecule in place into a small cluster of ~1,000 identical copies, all physically clustered within a few microns. This clonality is essential: SBS reads a cluster's average fluorescent signal each cycle, so you need many identical copies lighting up in sync to get a detectable signal above background.

Sequencing by synthesis. Each cycle, all clusters simultaneously incorporate one fluorescently labeled, reversibly-terminated nucleotide. A camera images the whole flowcell, recording which of four colors (or two, on two-channel chemistry used by NovaSeq/NextSeq) each cluster shows, then a chemical step removes the terminator and label so the next cycle can happen. Repeat for the length of the read (commonly 100-300 cycles). Because every cluster is imaged every cycle, you get massively parallel, simultaneous reads from millions to billions of clusters — this is "throughput" in the Illumina sense.

Paired-end sequencing. After reading from one end (Read 1), the library is reversed in place (not re-loaded) and the other end of the same molecule is read (Read 2). This gives you two reads per fragment, a known orientation (facing each other, "FR" orientation) and, because you know the library's size-selected fragment length distribution, an implied insert size — useful for structural variant detection and for resolving repeats, because you now have "anchor" information from both ends of a molecule even if one end lands in a repeat.

Index (barcode) reads. To run many samples in one lane, a short DNA barcode is ligated into the adapter next to the sample fragment. A separate short sequencing read (the index read, typically 6-10 bases, i7 and sometimes i5 for dual indexing) reads just the barcode so reads from different samples can be sorted out computationally (demultiplexing, section 4.3). Dual indexing (unique i7 AND i5 per sample, not just unique combinations) is the current best practice because it protects against index hopping (below).

Patterned vs. non-patterned flowcells. Non-patterned flowcells (older MiSeq/HiSeq models) have surface oligos randomly distributed; clusters form wherever a molecule happens to land, and cluster density is limited because clusters landing too close together overlap and become unresolvable ("polyclonal" or merged clusters). Patterned flowcells (NovaSeq, NextSeq 1000/2000, some iSeq) etch billions of nanowells at fixed, regular positions, one cluster per well by design. This raises achievable cluster density enormously (more data per flowcell, lower cost per Gb) but introduces a new failure mode: optical duplicates. Because wells are so close together, light from one well's cluster can bleed into the adjacent well's imaging window, and the base-calling software can occasionally split a single physical cluster's signal into two adjacent well calls, or two separate DNA molecules can seed adjacent wells during loading and be mistakenly treated downstream as duplicate reads from the same molecule. These are distinguished from PCR duplicates (library preparation amplified the same original fragment into >1 physical molecule before clustering) by their close spatial coordinates on the flowcell — a duplicate-marking tool that only looks at read coordinates (chromosome, position) cannot tell PCR duplicates from biologically-independent reads that happen to start at the same place, but tools aware of flowcell tile coordinates (e.g., Picard's MarkDuplicates with optical duplicate detection, or samblaster) can flag optical duplicates specifically by checking if the "duplicate" pair also sits physically adjacent on the flowcell.

Index hopping. On patterned flowcells using ExAmp chemistry (most current high-throughput Illumina instruments), free-floating un-ligated index primers or adapter fragments in the pooled library can, during the clustering chemistry, swap a cluster's index read to a different sample's barcode mid-process. The practical effect is that a small percentage (often 0.1-2%, higher with poor library cleanup) of reads get assigned to the wrong sample. This is why unique dual indexing (every sample gets a unique i7 AND a unique i5, not a reused combination) matters: with combinatorial (non-unique) dual indexing, a hopped read can still produce a valid-looking index combination belonging to another real sample in the pool and silently contaminate it; with unique dual indexing, a hopped read produces an index combination that does not match any expected sample and gets discarded during demultiplexing instead of silently contaminating a result.

4.1.3 Ion Torrent (semiconductor sequencing)

Ion Torrent also does sequencing by synthesis, but detects base incorporation by pH change, not fluorescence. Each of the four bases is flowed across the chip in turn (not all at once); when a base is complementary to the next unpaired template position, polymerase incorporates it and releases a hydrogen ion, which a semiconductor sensor under each well detects as a voltage change. Multiple identical bases in a row (a homopolymer) release a proportionally larger pH signal in one flow, and the chip must estimate "how many" from signal amplitude rather than counting discrete cycles. This estimation is the root of Ion Torrent's defining error mode: homopolymer-associated insertion/deletion errors — the chip over- or under-calls the length of homopolymer runs. Ion Torrent is fast and avoids the optical imaging bottleneck of Illumina, which made it historically attractive for same-day clinical and targeted panels, but its indel-heavy error profile in homopolymers is a persistent complication for variant calling.

4.1.4 PacBio HiFi (Circular Consensus Sequencing)

PacBio's long-read chemistry (single-molecule real-time, SMRT, sequencing) watches a single DNA polymerase, immobilized at the bottom of a tiny well (zero-mode waveguide), incorporate fluorescently labeled nucleotides in real time, one at a time, from a single template molecule — no clonal amplification, no clusters. Read length is long (often 10-25 kb) but raw single-pass accuracy is only moderate, with substitution and indel error scattered fairly uniformly.

HiFi (High-Fidelity) sequencing fixes the accuracy problem by ligating the template into a circular molecule (a SMRTbell) and letting the polymerase go around it repeatedly, generating multiple raw sub-reads of the same insert. Software then computes a consensus sequence across all passes of one molecule — this process is called Circular Consensus Sequencing (CCS). Because the passes are independent observations of the same template, random per-pass errors cancel out in the consensus, while any error that is systematic (rare, but exists for certain motifs) does not. The result: reads that are both long (10-20 kb typical insert size) and highly accurate (Phred-scored quality often >Q30, frequently described as ">99.9% accuracy" when enough passes are achieved). The practical cost is that circularizing and reading a molecule multiple times reduces total unique-molecule throughput compared to non-CCS long-read modes — you trade some raw throughput for consensus accuracy.

4.1.5 Oxford Nanopore (ONT)

Nanopore sequencing threads a single DNA (or RNA) strand through a protein nanopore embedded in a membrane, with a voltage applied across the membrane. As bases pass through the pore's narrowest constriction, they perturb the ionic current in a sequence-dependent way. The instrument records this current trace at very high temporal resolution — the "squiggle" — as a continuous, noisy analog signal, not discrete base calls. A separate computational step, basecalling, uses a trained neural network (historically recurrent networks, now mostly transformer-based models, e.g. ONT's Dorado basecaller) to translate the squiggle into a base sequence. Because basecalling is model-based and the model is regularly retrained and improved, the same raw squiggle data re-basecalled with a newer model can yield measurably higher accuracy — this is unusual among sequencing platforms and means archiving raw signal files (POD5/FAST5), not just FASTQ, has long-term value.

Nanopore offers two sequencing modes relevant to accuracy: simplex, reading one strand of the duplex DNA once, and duplex, where both strands of the same original double-stranded molecule are captured and sequenced back to back (because the two strands remain tethered through a hairpin adapter or are computationally paired), giving two independent observations of the same underlying sequence that can be consensus-combined — directly analogous in spirit to PacBio's CCS, but across the two strands of the duplex rather than multiple passes of a circle. Duplex reads reach accuracy comparable to HiFi; simplex reads are lower accuracy (historically substitution- and indel-prone, with indels concentrated in homopolymers and repetitive motifs, though per-base accuracy has improved substantially with newer pore chemistries and basecalling models). Read length on nanopore is essentially limited only by the integrity of the input DNA molecule — ultra-long protocols can produce reads of hundreds of kilobases to over a megabase, which is uniquely useful for resolving large structural variants, phasing haplotypes across long distances, and spanning highly repetitive regions no short-read or even standard long-read method can bridge.

4.1.6 Element Biosciences and Ultima Genomics, in one line each

Element (AVITI) uses a different clonal amplification and imaging chemistry than Illumina (avoiding some patented SBS steps) to produce short reads with an accuracy and cost profile broadly comparable to Illumina, marketed on lower reagent cost and instrument flexibility. Ultima Genomics uses a fundamentally different format — sequencing on an open, spinning silicon wafer rather than discrete flowcell lanes — to drive per-base cost down sharply for high-throughput short-read sequencing, at the cost of a somewhat different (currently less mature, improving) error profile and more specialized informatics tooling for wafer-level data.

4.1.7 Comparison table

A practical rule of thumb for choosing: if your question is "what variants does this genome/sample have, relative to a reference, at modest cost and scale," use Illumina (or Element/Ultima as cost-driven alternatives). If your question is "what is the actual structure of this genome, including repeats, large rearrangements, and haplotypes" or "what full-length transcript isoforms exist," use PacBio HiFi or Oxford Nanopore long reads, often both together (long reads for structure/phasing, short reads for cheap deep coverage of SNVs) in a hybrid design. If your question is "confirm this one specific base change in this one specific amplicon," Sanger is still the fastest, cheapest, most unambiguous answer.

4.2 Experimental design

Sequencing chemistry quality is wasted if the experiment answers the wrong question, or answers a question confounded with something you did not intend to measure. This section covers the design decisions that must be made before a single tube goes into a machine.

4.2.1 Biological vs. technical replication

A biological replicate is an independent biological unit: a different animal, a different patient, a different culture dish seeded and grown separately. A technical replicate is the same biological material measured more than once: splitting one extracted RNA sample into two library preps, or sequencing the same library on two lanes. Technical replicates estimate measurement noise (library prep variability, sequencing noise). Biological replicates estimate biological variability plus measurement noise, and only biological replicates let you generalize a claim ("treatment X changes expression of gene Y") beyond the single organism or sample you happened to measure. A common and serious design error is to run three technical replicates of one biological sample per condition and treat the resulting statistical significance as if it reflects biological reproducibility — it does not; it only shows the assay is precise, not that the biological effect is real or general. Differential expression and most other comparative statistical tests (Module 7, Statistics for Omics) require biological replicates as the unit of replication; technical replicates can be summed or averaged into one value per biological sample before that test, not treated as additional independent data points.

4.2.2 Depth and coverage mathematics

Coverage (also "depth") at a genomic position is the number of independent sequencing reads that overlap that position. Genome-wide average coverage is:

$$C = \frac{N \times L}{G}$$

where $N$ is the number of reads sequenced, $L$ is the read length (or total bases per read pair), and $G$ is the genome (or target) size in bases. This formula says: total bases sequenced, divided by the size of the thing you are covering, gives the average number of times each base is "seen."

Reads do not land uniformly — they land following, to a good first approximation, a Poisson process (random, independent starting positions). The Lander-Waterman model (originally derived for shotgun genome assembly) uses this to predict what fraction of the genome will be covered at least once, and the distribution of local depth, given an average coverage $C$. The probability that a specific base is covered by exactly $k$ reads is:

$$P(k) = \frac{C^k e^{-C}}{k!}$$

This is the Poisson probability mass function: $C$ is the mean of the distribution (the average coverage computed above), $k$ is a specific depth value you are asking about, and $e^{-C}$ normalizes the distribution so probabilities sum to 1. The shape says that even at a respectable average coverage, a meaningful fraction of bases will by chance be covered much less than average — the fraction with zero coverage is $P(0) = e^{-C}$. At $C = 10$, that is $e^{-10} \approx 0.0045\%$ of bases with no coverage at all — small but not zero, and this fraction grows fast as $C$ drops.

This directly drives practical depth requirements, because different downstream tasks need different minimum local, not just average, coverage:

A worked example: you want 30x average coverage of the human genome ($G \approx 3.1 \times 10^9$ bp, autosomes + sex chromosomes) using 150 bp paired-end reads (so $L = 300$ bp per read pair). Solve the coverage formula for $N$:

$$N = \frac{C \times G}{L} = \frac{30 \times 3.1 \times 10^9}{300} \approx 3.1 \times 10^8 \text{ read pairs}$$

That is about 310 million read pairs, or roughly 93 Gb of sequence (310M pairs × 300 bp). This is the number you hand to the sequencing core when requesting a run, or use to estimate cost.

4.2.3 Statistical power

Power is the probability of detecting a true effect of a given size, given your sample size, variability, and significance threshold. For RNA-seq differential expression, power depends on the number of biological replicates per group, the within-group variance (biological + technical noise), the fold change you want to detect, and the significance threshold corrected for the enormous number of simultaneous tests (one per gene — see Module 7 for multiple-testing correction). The practical consequence that surprises most new users: going from 3 to 6 biological replicates per group typically buys far more statistical power than doubling sequencing depth per sample, because biological variance, not sequencing noise, is usually the dominant source of uncertainty once depth is already adequate (the table above). If you have a fixed budget, spending it on more replicates at moderate depth almost always beats fewer replicates at very high depth.

4.2.4 Randomization and batch effects

A batch effect is systematic, non-biological variation introduced by when, where, or how samples were processed — different extraction days, different reagent lots, different technicians, different flowcells, different sequencing runs. Batch effects are dangerous specifically when they are confounded with the biological variable of interest: if all your "treatment" samples were processed on Monday and all your "control" samples on Tuesday, you cannot computationally separate "treatment effect" from "Monday-vs-Tuesday effect" — no normalization or statistical correction can recover information that was never collected. The fix is randomization: distribute samples from every biological condition across every batch, extraction day, library prep plate, and flowcell lane roughly evenly, so that batch variation is a nuisance variable you can model and remove (Module 7 covers batch-correction methods such as ComBat), not a variable perfectly aligned with your question.

A well-known real disaster illustrating this: a published microarray cancer study was later found to have sample processing dates nearly perfectly correlated with tumor subtype, so that the "molecular signature" the original analysis reported was substantially a batch artifact, not biology — this was a central example in later methodological critiques of batch effects in genomics (the "batch effect" literature following problems identified in some high-profile cancer genomics studies in the late 2000s, discussed at length by Leek and colleagues in their work on surrogate variable analysis). The general lesson generalizes directly to sequencing: always record and, where possible, randomize extraction batch, library prep batch, and flowcell/lane assignment against your biological groups before you generate data, not after.

4.2.5 Multiplexing, barcodes, and UMIs

Multiplexing means combining multiple samples' libraries, each carrying a distinct index barcode (section 4.1.2), into one sequencing run, then separating them computationally afterward (demultiplexing). This amortizes the fixed cost of a sequencing run across many samples and is standard practice for anything short of whole-genome-scale single-sample runs.

Unique Molecular Identifiers (UMIs) are short random sequences (commonly 6-12 bases) added to each original DNA/RNA molecule before PCR amplification, distinct from sample index barcodes. Every amplified copy of one original molecule carries the same UMI; different original molecules (almost certainly) get different random UMIs. This fixes a specific problem: PCR amplification is exponential and uneven, so some original molecules get massively over-amplified relative to others, and simply counting final reads overcounts the over-amplified ones. By grouping reads by UMI (plus mapping position) and collapsing each UMI group to a single consensus count or sequence, you recover an estimate of the number of original molecules, not the number of PCR copies — critical for accurate quantification in RNA-seq (especially single-cell) and for detecting rare mutations in deep sequencing (a true low-frequency variant must appear in multiple independent-UMI reads to be trusted; a variant appearing only within one UMI group is more likely a PCR or sequencing error introduced after the UMI was attached).

4.2.6 Controls

Positive controls are samples or spiked-in material with a known, expected result, included to confirm the entire pipeline (extraction through analysis) is working — e.g., a reference cell line with known variants, or a synthetic RNA spike-in (ERCC spike-ins) of known concentration added before library prep to check quantification linearity. Negative controls are samples expected to show no signal — a no-template control through library prep to catch reagent contamination, or a water/blank sample sequenced alongside real samples to catch index hopping or cross-contamination (any real-looking reads appearing in a water blank's demultiplexed output is a red flag). Omitting negative controls is how contamination and index-hopping problems go undetected until someone tries to explain an impossible result months later.

4.2.7 Avoiding confounded designs — a checklist

Before sequencing, ask: can I describe, for every sample, which batch/day/operator/reagent lot/flowcell lane it used, and is that assignment statistically independent of my biological variable of interest? If the honest answer is "no, treatment samples happen to all be in batch 1," the experiment cannot answer the question it was designed to answer, no matter how sophisticated the downstream computational correction — randomize the design before generating data, not after.

4.3 Raw data handling

4.3.1 From machine output to FASTQ: demultiplexing

Illumina sequencers write raw cluster intensity data as BCL files (one small binary file per cycle, per tile). bcl2fastq (legacy) and BCL Convert (current Illumina tool, bcl2fastq's replacement) convert these BCL files into per-sample FASTQ files, using the index reads and a sample sheet (a CSV mapping each expected i7/i5 index combination to a sample name) to sort reads by sample (demultiplexing) and discard or flag reads whose index doesn't match any expected combination within an allowed mismatch tolerance (typically 0 or 1 mismatch allowed, to avoid misassigning reads due to a single sequencing error in the index read itself).

# BCL Convert: demultiplex a NovaSeq run directory into per-sample FASTQs
bcl-convert \
  --bcl-input-directory /data/runs/240511_NovaSeq_Run42 \
  --output-directory /data/fastq/Run42 \
  --sample-sheet /data/runs/240511_NovaSeq_Run42/SampleSheet.csv \
  --bcl-sampleproject-subdirectories true
# Output: /data/fastq/Run42/<Sample_Name>/<Sample_Name>_S1_L001_R1_001.fastq.gz, _R2_001.fastq.gz
# plus Reports/ containing demultiplexing stats (reads per sample, % reads with
# unknown/unassigned index — this is where you check for index hopping damage)

A FASTQ file's four-line-per-read structure is: header line (@ + read ID, including flowcell/tile/x-y coordinates — the coordinates used for optical duplicate detection later), sequence line, + separator line, and a quality line of the same length, where each character encodes a Phred quality score $Q = -10 \log_{10}(p)$ ($p$ being the estimated probability that this specific base call is wrong; a Q30 base has a 1-in-1000 chance of being wrong, Q40 a 1-in-10,000 chance), offset by 33 in ASCII (the "Phred+33" encoding used by all current platforms).

4.3.2 FASTQ QC with FastQC, module by module

FastQC runs a fixed battery of checks on a FASTQ file and produces an HTML report with per-module pass/warn/fail flags. Knowing what each module actually plots — and which failures are expected artifacts vs. real problems — is the core skill.

fastqc sample_R1.fastq.gz sample_R2.fastq.gz -o qc/fastqc/ --threads 4
# produces sample_R1_fastqc.html, sample_R1_fastqc.zip (and R2 equivalents)

4.3.3 Adapter and quality trimming

Trimming removes two things: residual sequencing adapter (when the DNA fragment is shorter than the read length, sequencing "runs through" the fragment into the adapter) and low-quality bases (usually at read ends, where chemistry degrades).

# fastp: fast, combines adapter detection+trimming, quality filtering, and its own
# QC report in one tool (commonly preferred for its speed and sane defaults)
fastp \
  --in1 sample_R1.fastq.gz --in2 sample_R2.fastq.gz \
  --out1 sample_trimmed_R1.fastq.gz --out2 sample_trimmed_R2.fastq.gz \
  --detect_adapter_for_pe \
  --qualified_quality_phred 15 \
  --length_required 36 \
  --json fastp_report.json --html fastp_report.html \
  --thread 4
# key defaults: trims bases below Q15 from read ends, drops reads shorter than 36bp
# after trimming, auto-detects adapter by overlap analysis of R1/R2 (no adapter
# sequence needs to be specified for standard paired-end Illumina libraries)

# cutadapt: more manual control, useful when the adapter/UMI structure is non-standard
cutadapt \
  -a AGATCGGAAGAGC -A AGATCGGAAGAGC \
  -q 20,20 --minimum-length 36 \
  -o out_R1.fastq.gz -p out_R2.fastq.gz \
  sample_R1.fastq.gz sample_R2.fastq.gz
# -a/-A: 3' adapter sequence for R1/R2 (this is the standard Illumina TruSeq adapter);
# -q 20,20: quality-trim both ends of both reads to Q20

Trimmomatic remains common in legacy pipelines (especially published RNA-seq protocols) and uses a Java-based sliding-window approach (SLIDINGWINDOW:4:20 trims once the average quality in a 4-base window drops below Q20), but fastp has largely replaced it in new pipelines for speed and built-in reporting.

4.3.4 Contamination screening

FastQ Screen and Kraken2 both answer "what organisms/sequences are actually in this FASTQ file," but differently. FastQ Screen aligns a subsample of reads against a curated panel of reference genomes (human, mouse, common vectors, common lab contaminants like E. coli and mycoplasma, rRNA, adapter sequences) and reports the percentage of reads matching each — fast and interpretable because the panel is small and curated. Kraken2 classifies every read (or a subsample) against a much larger k-mer database spanning broad taxonomic breadth (e.g., all of RefSeq bacteria/archaea/viruses) using exact k-mer matching against a pre-built index, giving taxonomic composition down to species level — the right tool when you don't know what might be there (microbiome, suspected novel contamination) rather than checking against a short list of usual suspects.

# FastQ Screen: quick check against a defined contaminant panel
fastq_screen --aligner bowtie2 --conf fastq_screen.conf sample_R1.fastq.gz
# output: % reads mapping to each reference in the panel (human, mouse, rRNA, E.coli, etc.)

# Kraken2: broad taxonomic classification
kraken2 --db /refs/kraken2_standard_db \
  --paired sample_R1.fastq.gz sample_R2.fastq.gz \
  --report kraken2_report.txt --output kraken2_output.txt
# report columns: % reads, reads covered, reads assigned directly, rank code, taxID, name

4.3.5 MultiQC aggregation and the full annotated pipeline

MultiQC scans a directory for the output files of many tools (FastQC, fastp, Kraken2, alignment metrics, and more) and compiles them into one navigable HTML report with per-sample comparison tables and plots — essential once you have more than a handful of samples, because reading forty individual FastQC reports by hand does not scale.

#!/usr/bin/env bash
# raw_data_qc_pipeline.sh — demultiplex-to-clean-FASTQ pipeline, one sample pair
set -euo pipefail

SAMPLE=sampleA
R1=fastq/${SAMPLE}_R1.fastq.gz
R2=fastq/${SAMPLE}_R2.fastq.gz
OUT=qc_pipeline_out/${SAMPLE}
mkdir -p "${OUT}"/{fastqc_raw,fastp,fastqc_trimmed,screen,kraken2}

# 1. QC on raw reads
fastqc "$R1" "$R2" -o "${OUT}/fastqc_raw" --threads 4

# 2. adapter + quality trimming
fastp --in1 "$R1" --in2 "$R2" \
  --out1 "${OUT}/fastp/${SAMPLE}_trim_R1.fastq.gz" \
  --out2 "${OUT}/fastp/${SAMPLE}_trim_R2.fastq.gz" \
  --detect_adapter_for_pe --qualified_quality_phred 15 --length_required 36 \
  --json "${OUT}/fastp/${SAMPLE}.fastp.json" --html "${OUT}/fastp/${SAMPLE}.fastp.html" \
  --thread 4

# 3. QC on trimmed reads (confirm trimming worked: adapter content should now be ~0)
fastqc "${OUT}/fastp/${SAMPLE}_trim_R1.fastq.gz" "${OUT}/fastp/${SAMPLE}_trim_R2.fastq.gz" \
  -o "${OUT}/fastqc_trimmed" --threads 4

# 4. contamination screen (curated panel, fast)
fastq_screen --aligner bowtie2 --conf fastq_screen.conf \
  --outdir "${OUT}/screen" "${OUT}/fastp/${SAMPLE}_trim_R1.fastq.gz"

# 5. broad taxonomic check (only needed if screen is ambiguous or organism unknown)
kraken2 --db /refs/kraken2_standard_db \
  --paired "${OUT}/fastp/${SAMPLE}_trim_R1.fastq.gz" "${OUT}/fastp/${SAMPLE}_trim_R2.fastq.gz" \
  --report "${OUT}/kraken2/${SAMPLE}.kreport" --output "${OUT}/kraken2/${SAMPLE}.kraken" \
  --threads 4

# 6. aggregate everything into one report
multiqc "${OUT}" -o "${OUT}/multiqc"
# open ${OUT}/multiqc/multiqc_report.html — check per-sample tables for:
#   - %GC, total sequences, %duplication (fastqc_raw)
#   - % reads surviving filtering, adapter-trimmed % (fastp)
#   - adapter content now flat at ~0% across read position (fastqc_trimmed)
#   - % reads assigned to the expected organism vs. unexpected contaminants (screen/kraken2)

4.4 Alignment and post-processing

Once FASTQ files are clean, the next step is mapping reads to a reference genome and preparing the result for variant calling or quantification.

4.4.1 Choosing an aligner: bwa-mem2 vs. minimap2

bwa-mem2 is the standard short-read DNA aligner, a performance-optimized rewrite of BWA-MEM, built for accurate local alignment of reads 70 bp-1 Mb (in practice almost always used for short Illumina reads) against a reference genome, tolerant of a moderate rate of mismatches and short indels. minimap2 is a general-purpose aligner that handles short reads but is the standard choice for long reads (PacBio HiFi, Oxford Nanopore) and for genome-to-genome or splice-aware alignment, using a different underlying seeding/extension algorithm (minimizer-based seeding, well-suited to long sequences with higher raw error rates).

# bwa-mem2: short-read DNA alignment, with read-group annotation at alignment time
bwa-mem2 index reference.fa   # one-time index build (creates .0123, .amb, .ann, .bwt.2bit.64, .pac)

bwa-mem2 mem \
  -t 8 \
  -R '@RG\tID:sampleA_L001\tSM:sampleA\tPL:ILLUMINA\tLB:libA\tPU:Run42.L001' \
  reference.fa \
  sampleA_trim_R1.fastq.gz sampleA_trim_R2.fastq.gz \
  | samtools sort -@ 8 -o sampleA.sorted.bam -
# -t 8: use 8 threads
# -R: read group string, attached to every read in this invocation (see below for why this matters)
# piping directly into samtools sort avoids writing a large unsorted intermediate SAM/BAM to disk

# minimap2: long-read alignment (PacBio HiFi)
minimap2 -t 8 -ax map-hifi reference.fa sampleA.hifi.fastq.gz \
  | samtools sort -@ 8 -o sampleA.hifi.sorted.bam -
# -a: output SAM (vs. default PAF); -x map-hifi: preset tuned for HiFi error/length profile
# other common presets: map-ont (Nanopore simplex), map-pb (older PacBio CLR), splice (RNA)

Read groups (the -R string above) tag every read with metadata: ID (a unique identifier for this specific sequencing run/lane, used by duplicate marking and BQSR to know which reads came from the same physical sequencing event, since error profiles and PCR duplicates are run-specific), SM (sample name — the biological identity, used by variant callers to know which reads belong together for genotyping), PL (platform, e.g. ILLUMINA — used by some tools to select platform-appropriate error models), LB (library — distinguishes libraries if the same sample was prepped more than once, important for correctly identifying PCR duplicates, which can only occur within the same library), PU (platform unit, typically flowcell+lane — the most specific identifier of a single sequencing event). Getting read groups right matters enormously the moment you combine data from multiple runs, lanes, or libraries for the same sample (a routine situation): without correct ID/LB separation, duplicate marking and recalibration will either over-merge distinct runs (treating a cross-lane coincidental same-position pair as a PCR duplicate when it isn't) or under-merge them (missing real duplicates within a library that got split across lanes).

4.4.2 Sorting, indexing, duplicate marking

Aligned reads come out of the aligner in roughly the order reads were produced, not genome order. Most downstream tools (variant callers, samtools mpileup, visualizers like IGV) require coordinate-sorted BAM (Binary Alignment Map, the compressed binary form of SAM, Sequence Alignment/Map) files, and all of them require an index to do random access without reading the whole file.

# if not already sorted during alignment (see pipe above), sort now:
samtools sort -@ 8 -o sampleA.sorted.bam sampleA.bam
samtools index sampleA.sorted.bam        # produces sampleA.sorted.bam.bai

# mark duplicates (does not remove reads, flags them with SAM flag 0x400)
gatk MarkDuplicates \
  -I sampleA.sorted.bam \
  -O sampleA.dedup.bam \
  -M sampleA.dup_metrics.txt \
  --CREATE_INDEX true
# equivalent alternative: samtools markdup (faster, less memory, no per-library metrics file)

Duplicate marking flags reads that are PCR or optical duplicates (Module 4.1): read pairs whose outer coordinates match exactly are assumed to derive from one original DNA fragment amplified multiple times, not independent sampling of the genome. The tool does not delete these reads — it sets a SAM flag so downstream tools can ignore them — because duplicates are still useful for some QC (library complexity) and because deleting data silently is bad practice. Failing to mark duplicates before variant calling inflates apparent read depth and can produce false-positive variant calls supported only by repeated copies of one original molecule, not independent evidence. sampleA.dup_metrics.txt reports the duplication rate: under 5-10% is typical for a well-made WGS library at moderate input DNA; above 20-30% usually means low input DNA, over-amplification, or a degraded/low-complexity library, and should prompt a look at the library prep before trusting variant calls from that sample.