Agricultural Statistics & Data Science · Statistical genetics and biometrics · also Bioinformatics (sub-subject 7.3), Computational Biology
The third sub-subject of the ICAR Statistical Sciences group: molecular-biology basics, sequence databases, pairwise and multiple alignment, genome assembly, phylogenetics, and the statistics behind scoring and tree reliability. Taught from the statistician's side of the bench.
In the exams: AICE-JRF subject 69 Bioinformatics Units 1-5 (heavy: the whole paper); AIEEA-PG Group 07 sub-subject 7.3 (heavy); AIEEA-PG Group 07 bioinformatics databases (present)
Bioinformatics is the use of computation and statistics to store, retrieve, compare and interpret biological sequences and structures. In agriculture it underpins the sequencing of crop, livestock, pathogen and pest genomes, the discovery of genes and markers for breeding, the identification of pathogens and their resistance genes, and the study of how species and varieties are related. A statistician brings to it the scoring systems, probability models and resampling methods on which alignment, database search, assembly and tree-building depend.
This manual is written from the statistician's side of the bench. It begins with the molecular biology needed to understand what a sequence is and how genes are expressed and manipulated. It then covers the international sequence databases and file formats, pairwise alignment by dynamic programming, fast database search with BLAST and FASTA and the substitution matrices and statistics behind their scores, multiple sequence alignment, sequencing technologies and genome assembly, the measures used to judge an assembly and the coverage and error of reads, the construction of phylogenetic trees by parsimony, distance and likelihood methods, the assessment of tree reliability by bootstrap and jackknife, and an introduction to protein structure and functional annotation.
The subject is examined as a whole paper in agricultural PhD fellowship examinations in bioinformatics and as a sub-subject of postgraduate entrance examinations in the statistical sciences, where questions test definitions (database names, file formats, the genetic code), algorithms (filling a dynamic-programming matrix, the steps of BLAST, UPGMA and neighbour-joining), and simple computations (E-values, Phred scores, coverage, N50, Jukes-Cantor distances, numbers of possible trees). The manual names methods rather than software products, except where a database or tool name is itself the standard term.
Deoxyribonucleic acid (DNA) is a polymer of nucleotides, each made of a deoxyribose sugar, a phosphate and one of four bases: adenine (A), guanine (G), cytosine (C) and thymine (T). A and G are purines; C and T are pyrimidines. Watson and Crick's model of 1953 describes DNA as a right-handed double helix of two antiparallel strands held by hydrogen bonds between complementary bases: A pairs with T by two hydrogen bonds and G with C by three. Hence Chargaff's rules, A = T and G = C in double-stranded DNA, and the higher thermal stability of GC-rich DNA. Sequences are written 5′ to 3′, the direction in which polymerases add nucleotides. Replication is semi-conservative: each new double helix contains one old and one new strand, as Meselson and Stahl showed.
Ribonucleic acid (RNA) uses ribose and uracil (U) in place of thymine and is usually single-stranded. The central dogma describes the flow of information: DNA is replicated, transcribed into RNA, and messenger RNA (mRNA) is translated into protein. Reverse transcription (RNA to DNA), used by retroviruses and in the laboratory to make complementary DNA (cDNA), is an exception.
The genetic code reads mRNA in non-overlapping triplets (codons). Of the 64 codons, 61 specify the 20 amino acids and 3 (UAA, UAG, UGA) are stop codons; AUG codes for methionine and is the usual start codon. The code is degenerate (most amino acids have several codons, differing mainly in the third position) and nearly universal, with small variations in mitochondria and some microorganisms. Any DNA sequence has six possible reading frames, three on each strand; an open reading frame (ORF) is a stretch without stop codons beginning with a start codon, the starting point of gene finding.
GC content, the proportion of G and C bases, varies between and within genomes and affects melting temperature. For short oligonucleotides, the Wallace rule estimates the melting temperature as Tm ≈ 2(A + T) + 4(G + C) °C, used in primer design.
Transcription is carried out by RNA polymerase, which binds a promoter upstream of a gene and synthesises RNA complementary to the template strand. In bacteria, genes are often organised in operons transcribed together; the lac operon, described by Jacob and Monod, is the classic example of regulation by a repressor and an inducer. In eukaryotes, the primary transcript is processed: a 5′ cap is added, a poly(A) tail is added at the 3′ end, and introns are removed by splicing, leaving exons joined in the mature mRNA. Alternative splicing produces several proteins from one gene. Translation occurs on ribosomes, where transfer RNAs (tRNAs) carrying amino acids pair their anticodons with mRNA codons. Gene expression is regulated at every step, by transcription factors binding promoters and enhancers, by chromatin state and DNA methylation, by small RNAs, and by protein turnover.
Recombinant DNA technology joins DNA from different sources. Type II restriction enzymes cut DNA at specific, usually palindromic recognition sites; EcoRI, for example, cuts GAATTC between G and A, leaving single-stranded 'sticky' ends that pair with any fragment cut by the same enzyme. DNA ligase joins the fragments. Vectors carry inserted DNA into host cells: plasmids for small inserts, bacterial artificial chromosomes (BACs) for inserts of 100 kb or more. Transformation introduces the vector into a host, and selectable markers identify transformed cells. In plants, Agrobacterium-mediated transformation and particle bombardment introduce genes into crops.
The polymerase chain reaction (PCR) amplifies a target sequence between two primers through repeated cycles of denaturation (about 94-95 °C), primer annealing (typically 50-65 °C, depending on primer melting temperature) and extension (about 72 °C) by a heat-stable DNA polymerase. Each cycle can double the number of target molecules, so n cycles give up to 2ⁿ copies: 30 cycles give about 10⁹. Quantitative (real-time) PCR monitors amplification each cycle and estimates initial template amount from the threshold cycle, the basis of many expression and pathogen-detection assays. Gene editing with programmable nucleases, notably CRISPR-Cas systems, now allows targeted changes in crop and livestock genomes, and its regulation differs between countries.
Three databases together form the International Nucleotide Sequence Database Collaboration (INSDC): GenBank, maintained by the National Center for Biotechnology Information (NCBI) in the United States; the European Nucleotide Archive (ENA), maintained by the European Bioinformatics Institute (EMBL-EBI), which grew from the EMBL Data Library; and the DNA Data Bank of Japan (DDBJ). They exchange data daily, so a sequence submitted to any of them appears in all three with the same accession number. They are primary (archival) databases: they hold what submitters deposit, including redundant and sometimes erroneous entries.
Each entry has an accession number, a stable identifier, and a version number that increases when the sequence changes (for example an accession followed by .2). Derived, curated databases add value: RefSeq, at NCBI, provides a non-redundant reference set of genomic, transcript and protein sequences for each organism; raw sequencing reads go to the Sequence Read Archive and its European and Japanese counterparts; genome assemblies and their annotations are distributed through genome browsers such as Ensembl, which has divisions for plants, fungi, protists and metazoa, and through crop- and livestock-specific databases.
Protein sequences are held in UniProt, a collaboration of the EBI, the Swiss Institute of Bioinformatics (SIB) and the Protein Information Resource. UniProtKB has two sections: Swiss-Prot, manually reviewed and annotated, and TrEMBL, automatically annotated translations of nucleotide entries awaiting review. ExPASy, the SIB's portal, provides proteomics tools such as calculation of molecular weight and isoelectric point, translation, and access to PROSITE motifs. Protein structures are archived in the Protein Data Bank (PDB), managed by the worldwide PDB partnership. Domain and family databases (Pfam, now part of InterPro), pathway databases (KEGG and others) and the Gene Ontology support functional annotation.
NCBI's Entrez system links its databases (nucleotide, protein, gene, genome, structure, literature) so that a search in one leads to related records in others; the EBI and DDBJ offer equivalent search systems. Statisticians increasingly use programmatic access to download large data sets reproducibly rather than by manual web searches.
Three text formats are used constantly. The FASTA format consists of a header line beginning with '>' followed by an identifier and optional description, then one or more lines of sequence in single-letter codes. It is the input format of most alignment and search tools. Nucleotides use A, C, G, T (or U) and IUPAC ambiguity codes, such as N for any base, R for A or G (purine) and Y for C or T (pyrimidine); proteins use the 20 one-letter amino-acid codes.
The GenBank flat-file format (and the similar EMBL format) carries annotation: a LOCUS line with name, length, molecule type and date; DEFINITION, ACCESSION, VERSION and SOURCE/ORGANISM lines; a FEATURES table listing genes, coding sequences (CDS) with their coordinates and translations, and other features; and the sequence itself after ORIGIN. Coordinates are 1-based and inclusive, and 'complement(…)' marks features on the reverse strand.
The FASTQ format stores sequencing reads with base qualities, in four lines per read: a line beginning with '@' and the read identifier; the sequence; a line beginning with '+' (optionally repeating the identifier); and a quality string of the same length as the sequence, in which each character encodes the Phred quality of the corresponding base (character code minus an offset, usually 33).
Alignments of reads to a reference are stored in the SAM text format and its compressed binary equivalent BAM, with mapping positions, alignment details (CIGAR strings) and mapping quality; variants are stored in the Variant Call Format (VCF), with one line per variant and genotype columns per sample; genome annotations in GFF/GTF; and intervals in BED files, which, unlike GenBank, use 0-based start coordinates. Coordinate conventions (0- or 1-based, inclusive or half-open) are a frequent source of off-by-one errors when data move between formats.
An alignment arranges two sequences in columns so that residues thought to derive from a common ancestor (homologous positions) are aligned, inserting gaps to represent insertions or deletions (indels). Its quality is measured by a score: the sum of substitution scores for aligned residue pairs minus penalties for gaps. For DNA, a simple scheme scores matches positively and mismatches negatively (for example +1 and −1), and transitions (purine to purine or pyrimidine to pyrimidine) may be penalised less than transversions. For proteins, substitution matrices such as PAM and BLOSUM give log-odds scores for each pair of amino acids (Chapter 4).
Gap penalties may be linear, a cost d per gap position (a gap of length g costs gd), or affine, a larger cost for opening a gap and a smaller cost for each extension, g(d) = −(o + e(g − 1)) or similar. Affine penalties reflect the biology: one indel event often removes or inserts several residues together, so one long gap is more plausible than several short ones.
Two descriptive measures summarise an alignment: identity, the percentage of aligned positions with identical residues, and similarity, the percentage with identical or conservatively substituted residues. Identity depends on how gaps and alignment length are counted, so it should be reported with its definition. Homology is a yes-or-no statement about common ancestry; sequences are homologous or not, never 'percent homologous'. Homologues that diverged by speciation are orthologues; those that diverged by gene duplication are paralogues.
A dot plot, a matrix with one sequence along each axis and a dot where residues (or short windows) match, gives a visual overview: diagonals show similar regions, parallel diagonals show repeats, and broken diagonals show indels. It is a useful first look but does not give an optimal alignment.
The number of possible alignments of two sequences grows exponentially with their length, so they cannot all be scored. Dynamic programming finds the optimal alignment exactly in time proportional to the product of the lengths, by building a matrix F in which F(i, j) is the best score of an alignment of the first i residues of sequence x with the first j residues of sequence y.
The Needleman-Wunsch algorithm (1970) finds the optimal global alignment, covering both sequences from end to end. With a linear gap penalty d, the matrix is initialised with F(0, 0) = 0, F(i, 0) = −id and F(0, j) = −jd, and filled by the recurrence F(i, j) = max{F(i − 1, j − 1) + s(xᵢ, yⱼ), F(i − 1, j) − d, F(i, j − 1) − d}, where the three terms correspond to aligning xᵢ with yⱼ, aligning xᵢ with a gap, and aligning yⱼ with a gap. The optimal score is F(m, n), in the bottom-right cell, and the alignment is recovered by a traceback, following the pointers that recorded which term gave each maximum.
The Smith-Waterman algorithm (1981) finds the optimal local alignment, the best-scoring pair of segments, which is what is wanted when sequences share only a domain or when a short sequence is sought within a long one. It modifies the recurrence by adding a fourth option, zero: F(i, j) = max{0, F(i − 1, j − 1) + s(xᵢ, yⱼ), F(i − 1, j) − d, F(i, j − 1) − d}, initialises the first row and column to zero, and starts the traceback at the highest value anywhere in the matrix, ending when a zero is reached. For local alignment to be meaningful, the expected score of aligning random residues must be negative; otherwise long random alignments would score highly.
Both algorithms need time and memory proportional to mn. For two genes of a few thousand bases this is trivial, but for searching a whole database it is too slow, which led to the heuristic methods of the next chapter. Variants include semi-global (overlap) alignment, which does not penalise end gaps and is used to detect overlaps between sequencing reads, and affine-gap versions (Gotoh's algorithm) that keep three matrices.
A substitution matrix gives a score for aligning each pair of residues. Modern matrices are log-odds matrices: the score for residues a and b is s(a, b) = (1/λ) log[q_ab/(p_a p_b)], where q_ab is the frequency with which a and b are aligned in true alignments of related proteins, p_a and p_b are background frequencies, and λ is a scaling constant (scores are rounded to integers). A positive score means the pair occurs in related proteins more often than by chance; a negative score, less often.
PAM (point accepted mutation) matrices were derived by Margaret Dayhoff and colleagues in 1978 from alignments of closely related proteins. One PAM unit corresponds to 1% of amino acids changed (an accepted point mutation per 100 residues). The PAM1 transition matrix was estimated from closely related sequences and extrapolated to greater distances by multiplying it by itself: PAM250 = PAM1²⁵⁰, which corresponds to sequences of roughly 20% identity. Higher PAM numbers suit more distant sequences. The extrapolation assumes the substitution process is the same at all distances, which is a weakness.
BLOSUM (blocks substitution matrix) matrices were derived by Henikoff and Henikoff in 1992 directly from conserved, ungapped blocks of aligned protein families. For BLOSUM62, sequences in each block that were at least 62% identical were clustered and down-weighted, so that the matrix reflects substitutions between sequences less than 62% identical. Higher BLOSUM numbers suit closer sequences, the reverse of PAM: BLOSUM80 for close relatives, BLOSUM45 for distant ones. BLOSUM62 is the default for most protein database searches. Approximate correspondences often quoted are PAM250 with BLOSUM45 and PAM120 with BLOSUM80.
DNA scoring usually uses simple match and mismatch scores, sometimes distinguishing transitions from transversions.
Database searching compares a query sequence with millions of entries, so speed matters more than guaranteed optimality. FASTA (Lipman and Pearson, 1985, and Pearson and Lipman, 1988) first finds short exact matches of length ktup (word length, often 1-2 for proteins and 4-6 for DNA) between query and database sequence, identifies diagonals with many such matches, joins nearby high-scoring regions, and finally performs a banded Smith-Waterman alignment around the best region. The FASTA format of sequence files takes its name from this program.
BLAST (Basic Local Alignment Search Tool; Altschul and colleagues, 1990) is the most used search method. It (1) breaks the query into short words (default length 3 for proteins and 11 for standard nucleotide searches, longer for fast nucleotide modes) and lists all words that score at least a threshold T against each query word using the substitution matrix (the neighbourhood); (2) scans the database for exact hits to these words (seeds); (3) extends each hit in both directions without gaps, and later versions with gaps, keeping it while the score stays within a drop-off value of its maximum; and (4) reports high-scoring segment pairs (HSPs) with their statistics. The two-hit method, requiring two nearby hits on the same diagonal before extension, speeds up later versions.
BLAST programs differ by query and database type: blastn (nucleotide against nucleotide), blastp (protein against protein), blastx (translated nucleotide query against proteins, useful for finding the protein coded by a new DNA sequence), tblastn (protein query against translated nucleotide database, useful for finding unannotated genes in genomes) and tblastx (translated against translated). Position-specific iterated BLAST builds a profile from the first round's hits and searches again, detecting more distant relatives.
Low-complexity regions (such as runs of one amino acid) and repeats produce many spurious hits and are masked by default filters.
The statistics of local alignment scores were worked out by Karlin and Altschul (1990). For ungapped local alignments of random sequences, the number of HSPs with score at least S expected by chance is E = K m n e^(−λS), where m and n are the lengths of the query and the database (total residues), and K and λ are constants that depend on the scoring system and residue composition. Gapped alignments follow the same form approximately, with K and λ estimated by simulation. The scores of the best alignments follow an extreme-value (Gumbel) distribution, not a normal one.
Because raw scores depend on the scoring system, BLAST reports a normalised bit score, S′ = (λS − ln K)/ln 2, in terms of which E = m n 2^(−S′). The E-value is the expected number of chance hits with at least this score in a database of this size. An E-value of 10⁻⁵⁰ means essentially no chance of such a hit by chance; an E-value of 1 means one such hit is expected by chance alone. E-values scale with database size: the same alignment against a database ten times larger has an E-value ten times larger. For small E, the probability of at least one chance hit is P = 1 − e^(−E) ≈ E.
Example: a query of 250 residues searched against a database of 10⁸ residues gives a hit with bit score 50. Then E = 250 × 10⁸ × 2⁻⁵⁰ ≈ 2.5 × 10¹⁰ × 8.9 × 10⁻¹⁶ ≈ 2.2 × 10⁻⁵, a significant hit. For the same bit score against a database 100 times larger, E ≈ 2.2 × 10⁻³.
Common practice treats E-values below about 10⁻³ to 10⁻⁵ as evidence of homology for proteins, though thresholds depend on purpose. A significant E-value supports homology; a non-significant one does not prove absence of homology, since distant homologues may score poorly.
A multiple sequence alignment (MSA) aligns three or more homologous sequences so that each column contains residues derived from a common ancestral position. MSAs are the input to phylogenetic analysis, the source of conserved regions for primer and probe design, the basis of profiles and hidden Markov models for detecting distant family members, and evidence for functional and structural constraints (conserved columns often mark active sites or structural cores).
The natural extension of dynamic programming to N sequences of length L requires time and memory proportional to Lᴺ, which is feasible only for a handful of short sequences. The usual objective, the sum-of-pairs (SP) score, adds the pairwise scores of all pairs of sequences in each column. Because exact optimisation is impractical, heuristic methods are used.
Progressive alignment, the most widely used approach (as in the long-established Clustal family of programs), has three steps: compute pairwise distances between all sequences (from pairwise alignments or fast k-mer counts); build a guide tree from the distances, by neighbour-joining or UPGMA; and align sequences progressively following the tree, starting with the most similar pair and adding sequences or groups (profiles) in order of decreasing similarity. Its weakness is that errors made early are never corrected ('once a gap, always a gap'), and the result depends on the guide tree.
Iterative refinement methods realign subsets of sequences repeatedly and accept changes that improve the objective score. Consistency-based methods build a library of pairwise alignments and favour placements that agree across many pairs, which improves accuracy on difficult sets at a computational cost. Fast methods based on k-mer distances and fast Fourier transforms scale to tens of thousands of sequences. Accuracy is benchmarked against reference alignments based on known protein structures.
An MSA is summarised by its consensus sequence, the most common residue in each column, and more fully by a profile or position-specific scoring matrix (PSSM), the frequency (or log-odds score) of each residue in each column. A sequence logo displays the information content of each column as a stack of letters whose heights are proportional to their frequencies and whose total height is the column's information content in bits, log₂ 4 − H for DNA (maximum 2 bits) or log₂ 20 − H for proteins (maximum about 4.32 bits), where H is the Shannon entropy of the column.
Profile hidden Markov models (HMMs) represent a family alignment as a chain of match, insert and delete states with position-specific emission and transition probabilities. They are trained from an alignment and used to search databases for further members and to align new sequences to the family. Protein domain databases, such as Pfam, are collections of profile HMMs built from curated seed alignments.
Alignments must be checked before use. Poorly aligned regions, often at the ends or in variable loops, contain columns that are not truly homologous and can mislead phylogenetic inference; such columns are often trimmed, by eye or by automatic filters that remove columns with many gaps or low consistency. Protein-coding DNA is better aligned at the protein level and then back-translated, so that gaps fall between codons. Alignment uncertainty can be propagated into downstream analyses by using several alternative alignments or by methods that co-estimate alignment and tree, but most analyses still take one alignment as fixed, a source of overconfidence that a careful analyst should acknowledge.
Sanger sequencing, the chain-termination method developed by Frederick Sanger in 1977, synthesises DNA copies of a template in the presence of dideoxynucleotides (ddNTPs), which lack the 3′ hydroxyl group and stop chain extension. Fragments ending at each position are separated by size, now by capillary electrophoresis with fluorescently labelled ddNTPs, and read in order. Sanger reads are about 700-1,000 bases long and highly accurate, and the method remains the standard for verifying individual sequences, but its throughput is low and its cost per base high.
Second-generation (short-read, massively parallel) sequencing, from the mid-2000s, sequences millions to billions of clonally amplified DNA fragments simultaneously. The dominant approach, sequencing by synthesis with reversible fluorescent terminators, reads one base per cycle on clusters of identical molecules, producing reads typically of 100-300 bases with low error rates, mostly substitutions, which rise towards the ends of reads. Paired-end sequencing reads both ends of each fragment, giving two reads at a known approximate distance, which helps mapping and assembly. The cost per base fell by orders of magnitude, making whole-genome resequencing of crop panels and livestock populations routine.
Third-generation (long-read, single-molecule) technologies sequence individual molecules without amplification. Single-molecule real-time sequencing observes a polymerase incorporating labelled nucleotides; nanopore sequencing measures changes in ionic current as DNA passes through a protein pore. Reads are typically tens of kilobases and can exceed a hundred kilobases. Raw error rates were high, mostly insertions and deletions, but circular consensus methods and improved chemistry and base-calling have brought accuracy close to short reads. Long reads span repeats that short reads cannot, transforming the assembly of large, repetitive plant genomes.
Choice of technology depends on purpose: short reads for resequencing, genotyping and expression at low cost; long reads for de novo assembly, structural variants and full-length transcripts; and combinations of both. Complementary scaffolding technologies, such as chromatin conformation capture (Hi-C) and optical maps, order and orient assembled sequences into chromosome-scale assemblies.
Genome assembly reconstructs a genome from reads much shorter than it. De novo assembly uses only the reads; comparative (reference-guided) assembly maps reads to the genome of a related individual or species and builds the new sequence from that alignment, which is faster but can miss sequence absent from the reference and inherit its errors.
The overlap-layout-consensus (OLC) approach computes pairwise overlaps between reads, builds an overlap graph whose nodes are reads and edges are overlaps, finds a layout (a path through the graph), and derives a consensus sequence. It suits long reads, for which the number of reads is manageable and overlaps are long. The greedy approach, which repeatedly merges the pair of reads with the best overlap, was used in early assemblers and is simple but easily misled by repeats.
The de Bruijn graph approach, which dominates short-read assembly, breaks every read into overlapping k-mers (substrings of length k); nodes represent (k − 1)-mers and edges represent k-mers, so that the genome corresponds to a path traversing every edge (an Eulerian path). It avoids computing all pairwise overlaps among billions of reads. The choice of k trades sensitivity (smaller k connects more reads but collapses repeats) against specificity (larger k resolves more repeats but needs higher coverage and is more affected by errors).
Repeats are the central difficulty in all approaches: a repeat longer than the reads or the k-mer size creates ambiguous paths, so assemblies break into contigs (contiguous sequences) at repeats. Paired-end and mate-pair reads, long reads, Hi-C and genetic or optical maps link contigs into scaffolds, with gaps of estimated size filled with Ns. Plant genomes are especially challenging: many are large (bread wheat about 16 gigabases), rich in transposable elements, and polyploid, with similar homoeologous subgenomes. Heterozygous genomes produce two haplotypes that assemblers may separate or collapse, and haplotype-resolved assembly is an active area. Genome size can be estimated before assembly from the k-mer frequency spectrum: total number of k-mers divided by the depth at the main peak.
Each base call carries a quality score on the Phred scale: Q = −10 log₁₀ P, where P is the estimated probability that the base call is wrong. Q10 means a 1 in 10 error probability (90% accuracy), Q20 means 1 in 100 (99%), Q30 means 1 in 1,000 (99.9%), and Q40 1 in 10,000. Reads are trimmed or filtered by quality before analysis, and adapters are removed. The average error rate of a set of reads is computed by converting each Q back to P and averaging, not by averaging Q values.
Coverage (depth) is the average number of reads covering each base: C = LN/G, where L is the read length, N the number of reads and G the genome size. Sequencing 40 million paired reads of 150 bases each (80 million reads in all) of a 400-megabase genome gives C = 150 × 80 × 10⁶/(400 × 10⁶) = 30×.
The Lander-Waterman model (1988) treats read start positions as random (Poisson) along the genome. The depth at a given base is then Poisson with mean C, so the probability that a base is not covered by any read is e^(−C): at 1× about 37% of the genome is missed; at 8×, about 0.03%. Real coverage is less uniform than Poisson because of GC bias, repeats and library preparation, so higher nominal coverage is needed. Typical targets are about 30× for calling variants in a diploid with short reads (enough to see both alleles of heterozygotes reliably) and higher depth or long reads for de novo assembly; requirements vary with technology, ploidy and purpose.
Duplicate reads from PCR amplification of the same original fragment are marked and usually removed, because they give false confidence in errors that arose during amplification. Variant callers combine base qualities, mapping qualities and depth in likelihood models of the genotypes; low depth, low mapping quality (reads that map equally well to several places) and strand bias are standard filters.
An assembly is judged on contiguity, completeness and correctness.
Contiguity is summarised by N50: sort contigs (or scaffolds) from longest to shortest and add their lengths until the running total reaches half of the total assembly length; the length of the contig at which this happens is the N50. L50 is the number of contigs needed to reach it. For contigs of 100, 80, 60, 50, 40, 30, 20, 10, 5 and 5 kb (total 400 kb), the cumulative sums are 100, 180, 240, …, so half (200 kb) is reached at the third contig and N50 = 60 kb, L50 = 3. NG50 uses half of the estimated genome size rather than the assembly size, which avoids rewarding assemblies that omit difficult regions. A higher N50 is better, but N50 can be inflated by joining contigs wrongly.
Completeness is assessed by the proportion of the estimated genome size that is assembled, by the fraction of reads that map back to the assembly, and by the presence of a benchmark set of genes expected to occur once in nearly all species of a lineage (conserved single-copy orthologues): an assembly that contains 98% of them complete and in single copy is likely nearly complete in gene space, while many duplicated copies suggest uncollapsed haplotypes or contamination. The proportion of known transcripts that align full-length is another check.
Correctness is assessed by mapping reads and looking for inconsistent coverage or pairs, by comparison with genetic maps or with independent long-range data, by k-mer comparisons between reads and assembly (which estimate base-level quality as a consensus quality value on the Phred scale), and by checking for contamination from other organisms (bacteria, fungi, the human sequencer) by searching against databases. Misassemblies, such as joins of sequences that are distant in the genome or collapsed repeats, are the most damaging errors because they mislead downstream mapping and gene annotation.
Assemblies are reported with their sequencing data, coverage, assembly method and versions, contiguity statistics, completeness scores and known limitations, and deposited in the public databases.
A phylogenetic tree represents the evolutionary relationships among sequences or taxa. Its external nodes (leaves) are the observed sequences; internal nodes are inferred ancestors; branches may have lengths proportional to the amount of change. A rooted tree has a node representing the common ancestor of all leaves and gives the direction of time; an unrooted tree shows only relationships. Trees are rooted by an outgroup (a sequence known to lie outside the group of interest) or by assuming a molecular clock. A clade is an ancestor with all its descendants.
The number of possible trees grows explosively. For n taxa, the number of unrooted bifurcating trees is (2n − 5)!! = 3 × 5 × … × (2n − 5), and of rooted trees (2n − 3)!!. Five taxa have 15 unrooted and 105 rooted trees; ten taxa have over two million unrooted trees, and fifty taxa more than there are atoms in the observable universe. Exhaustive search is possible only for about a dozen taxa; beyond that, heuristic searches rearrange trees (nearest-neighbour interchange, subtree pruning and regrafting) to climb towards the best.
Distance methods first compute a matrix of pairwise evolutionary distances and then build a tree from it. The observed proportion of differing sites, the p-distance, underestimates the true number of substitutions because of multiple hits at the same site. Model-based corrections estimate it: under the Jukes-Cantor model, in which all substitutions are equally likely, d = −(3/4) ln(1 − 4p/3). For p = 0.3, d = 0.383 substitutions per site. The Kimura two-parameter model allows different rates for transitions and transversions.
UPGMA (unweighted pair group method with arithmetic mean) repeatedly joins the two closest clusters and sets the distance from the new cluster to each other cluster as the average of the distances of its members. It produces a rooted, ultrametric tree, which is correct only if a molecular clock holds (equal rates in all lineages); otherwise it can give the wrong topology. Neighbour-joining (Saitou and Nei, 1987) does not assume a clock: at each step it joins the pair that minimises the total tree length, using a corrected criterion Q(i, j) = (n − 2)d(i, j) − Σₖd(i, k) − Σₖd(j, k), and gives an unrooted tree with unequal branch lengths. It is fast and widely used, and is consistent when distances are additive.
Maximum parsimony chooses the tree that requires the fewest character changes (substitutions) to explain the observed sequences. For each candidate tree, the minimum number of changes at each site is computed (by Fitch's algorithm for unordered characters), and the tree with the smallest total is preferred. Only parsimony-informative sites matter: a site is informative if at least two character states each occur in at least two taxa; sites that are constant, or in which only one taxon differs (singletons), give the same count on every tree. Parsimony is intuitive and model-free, but it can be statistically inconsistent: when two long branches are separated by a short one, parsimony tends to group the long branches together because of parallel changes, a bias called long-branch attraction.
Maximum likelihood (ML), developed for phylogenetics by Felsenstein (1981), computes the probability of the observed alignment given a tree (topology and branch lengths) and an explicit model of substitution, and chooses the tree and parameters that maximise it. The likelihood of each site is summed over all possible states at internal nodes efficiently by Felsenstein's pruning algorithm; sites are assumed independent, so the log-likelihood of the alignment is the sum over sites. Substitution models range from Jukes-Cantor (JC69; equal base frequencies and rates) through Kimura (K80; transition/transversion ratio), HKY (adding unequal base frequencies) to the general time-reversible model (GTR; six exchange rates and four frequencies), usually with rate variation among sites modelled by a gamma distribution (+Γ) and a proportion of invariable sites (+I). Models are compared by likelihood-ratio tests for nested models or by information criteria (AIC, BIC). ML is statistically consistent when the model is adequate, and is less affected by long-branch attraction than parsimony, but it is computationally demanding.
Bayesian phylogenetics uses the same likelihoods with prior distributions on trees and parameters and samples the posterior distribution by Markov chain Monte Carlo, giving posterior probabilities of clades directly. Molecular-clock dating with calibrations from fossils or known divergence events estimates divergence times, for example of domesticated crops from wild relatives.
A tree inferred from one alignment is an estimate, and its clades may be supported strongly or weakly by the data. Felsenstein introduced the non-parametric bootstrap to phylogenetics in 1985. A bootstrap replicate is made by sampling alignment columns (sites) with replacement until a new alignment of the same length is formed; some sites appear several times and others not at all. A tree is inferred from each replicate by the same method used for the original, and this is repeated many times (commonly 100 to 1,000). The bootstrap support of a clade is the percentage of replicate trees that contain it. Supports are displayed on the branches of the original tree or summarised in a majority-rule consensus tree, which contains every clade found in more than 50% of the replicates.
Interpretation needs care. Bootstrap support measures how consistently the data support a clade under the chosen method and model, a measure of repeatability; it is not the probability that the clade is true. A systematic error, such as long-branch attraction or a misspecified model, can give high bootstrap support to a wrong clade, because every replicate inherits the same bias. Conversely, simulation studies (notably by Hillis and Bull) suggested that bootstrap proportions are conservative under favourable conditions, with values around 70% or more often corresponding to a high probability of the clade being correct. Conventions differ, but values of 70% and above are commonly described as moderate to good support and values above 90-95% as strong; values below 50% are not shown in majority-rule trees.
The bootstrap assumes sites are independent and identically distributed, which is not strictly true (codon positions, structural constraints), and it is computationally expensive for ML, so fast approximations (such as ultrafast bootstrap and branch tests based on approximate likelihood ratios) are widely used; their values are not on the same scale as standard bootstrap percentages and must be labelled correctly.
The jackknife resamples by deletion rather than replacement: each replicate removes a fixed fraction of sites (often 50%, or about 37% to mimic the bootstrap, since a bootstrap sample omits on average a fraction 1/e ≈ 0.368 of sites) and rebuilds the tree from the remainder. Jackknife support is the percentage of replicates containing a clade and behaves similarly to the bootstrap; it is common in parsimony analyses. The delete-one jackknife, removing one site or one taxon at a time, is used to assess the influence of individual sites or taxa.
Consensus trees summarise sets of trees, from resampling or from several equally good trees. The strict consensus contains only clades present in all trees and is often poorly resolved. The majority-rule consensus contains clades present in more than half of the trees, with their frequencies; extended majority-rule adds compatible clades with lower frequencies to resolve the tree further.
Other measures include Bremer support (decay index) in parsimony, the number of extra steps needed before a clade is lost from the most parsimonious trees; Bayesian posterior probabilities, which tend to be higher than bootstrap proportions for the same clade and are interpreted as the probability of the clade given the data, model and priors; and statistical tests comparing whole topologies (for example the Kishino-Hasegawa, Shimodaira-Hasegawa and approximately unbiased tests), which ask whether an alternative tree is significantly worse than the best one, with corrections for selection bias.
In practice a careful analysis reports the method, model, number of replicates and the type of support value, shows support on every branch, and discusses weakly supported or conflicting clades rather than presenting the tree as certain. Different genes may support different trees because of incomplete lineage sorting, hybridisation or horizontal transfer; multi-gene and genome-scale analyses combine genes by concatenation or by coalescent-based species-tree methods.
Proteins are chains of amino acids joined by peptide bonds. Their structure is described at four levels. The primary structure is the amino-acid sequence. The secondary structure consists of local regular patterns stabilised by hydrogen bonds between backbone atoms: the α-helix, with 3.6 residues per turn and a rise of 1.5 Å per residue (5.4 Å per turn), and the β-sheet, formed by extended strands lying side by side in parallel or antiparallel orientation; turns and loops connect them. The tertiary structure is the three-dimensional fold of one chain, stabilised by hydrophobic interactions, hydrogen bonds, ionic interactions and disulphide bridges. The quaternary structure is the arrangement of several chains (subunits) in a complex.
Backbone conformation is described by two dihedral angles per residue, φ and ψ; the Ramachandran plot of ψ against φ shows the allowed regions, with clusters for α-helix and β-sheet, and is used to check the quality of experimental and predicted structures. Anfinsen's experiments showed that the amino-acid sequence contains the information needed to determine the fold, which is the basis of structure prediction.
Experimental structures are determined mainly by X-ray crystallography, nuclear magnetic resonance (NMR) spectroscopy and cryo-electron microscopy, and deposited in the Protein Data Bank. Prediction methods include homology (comparative) modelling, which builds a model from the structure of a homologous template and is reliable when sequence identity is high (accuracy falls sharply below about 30% identity); fold recognition (threading), which fits a sequence to known folds when similarity is too weak to detect by sequence; and ab initio methods. Since about 2020, deep-learning methods trained on known structures and multiple sequence alignments have predicted many structures with accuracy close to experiment, an advance recognised by the 2024 Nobel Prize in Chemistry, and databases of predicted structures now cover most known proteins. Predicted models carry per-residue confidence scores that must be considered. Structures are compared by superposition and the root-mean-square deviation (RMSD) of corresponding atoms.
Annotation assigns biological meaning to sequences: where the genes are, what they encode, and what the products do.
Structural annotation finds genes. Ab initio gene finders use statistical models, typically hidden Markov models, of coding sequence, splice sites, start and stop codons trained on known genes of the species; evidence-based methods align transcripts (RNA-seq, full-length cDNA) and proteins from related species to the genome; combined pipelines integrate both. Repeats and transposable elements are identified and masked first, because they are abundant in plant genomes and resemble genes. The quality of gene models varies widely between genome projects.
Functional annotation infers what a gene product does, mostly by homology: a significant BLAST hit to a characterised protein, conserved domains detected by profile HMMs (Pfam and other InterPro member databases), motifs and patterns (PROSITE), signal peptides and transmembrane helices predicted by dedicated models, and orthology to genes of model organisms. Functions are recorded with controlled vocabularies, chiefly the Gene Ontology, which has three independent ontologies: molecular function (what a product does at the molecular level), biological process (the larger process it contributes to) and cellular component (where it acts). Each annotation carries an evidence code showing whether it comes from experiment or from computational inference. Genes are also mapped to metabolic and signalling pathways.
Annotation by homology has limits that a statistician should keep in mind. Transfer of function from a best hit is only as good as the annotation of that hit, and errors propagate through databases as automatically annotated sequences are themselves used as references. Paralogues may have diverged in function, so orthology is a better basis than similarity alone. Many genes, even in well-studied crops, remain annotated only as 'hypothetical protein' or 'protein of unknown function'. Enrichment analyses, testing whether a list of genes (for example differentially expressed or GWAS candidates) contains more genes of a given GO term than expected, must use the correct background set and correct for multiple testing, because thousands of terms are tested.
A student preparing for a bioinformatics fellowship examination must fill a small Needleman-Wunsch or Smith-Waterman matrix by hand, choose the right BLAST program, compute E-values, Phred error rates, coverage, N50 and Jukes-Cantor distances, count possible trees, and run UPGMA and neighbour-joining on a small distance matrix. In thesis work the same student downloads sequences, aligns them, builds a tree with bootstrap support, and annotates candidate genes from a QTL or GWAS interval, reporting every database version and parameter used so that others can repeat it.
A bioinformatician five years into work at a research institute, genomics facility or breeding company runs sequencing projects end to end: chooses technology and depth, checks read quality, assembles or maps reads, calls variants, assesses assemblies with contiguity and completeness measures, annotates genes, and builds phylogenies of pathogens or germplasm. Such a person writes reproducible pipelines, documents versions and parameters, and explains to biologists what an E-value or bootstrap value does and does not mean, and when a result needs laboratory confirmation.
A senior bioinformatician with twenty-five years of experience leads a genomics platform or a national crop or livestock genome project: plans sequencing strategies for large and polyploid genomes, sets data standards and submission policies to the international databases, manages computing infrastructure and data sharing agreements, and judges when new technologies are mature enough to adopt. Such a person also coordinates with plant protection and animal health authorities on pathogen sequencing for surveillance and outbreak investigation in crops and livestock, and mentors staff.
A professor or research leader with a PhD and long experience develops methods for alignment, assembly, phylogenetics and annotation of complex agricultural genomes, including pangenomes that capture variation across many varieties, and integrates statistical modelling with sequence analysis. Such a person trains bioinformaticians, sets curricula linking statistics, computing and biology, advises international germplasm and genomics consortia, and promotes open data and reproducible analysis, including fair sharing of sequence data from the countries where crop diversity and its wild relatives originate.
The manual explains algorithms and statistics but does not teach any software's command syntax. Detailed models of sequencing error, variant calling algorithms, transcriptome assembly, metagenomics, pangenome graphs, structural-variant detection and molecular dynamics are not covered. Default parameters of search tools (word sizes, thresholds) and typical read lengths and accuracies change as technologies and software versions change; values given are typical at the time of writing and should be checked against current documentation. Conventions for interpreting bootstrap support and E-value thresholds differ between fields and are given as common practice, not rules.
Original practice questions written to the syllabus and pattern of the exams named in each question's examTags. They are not past papers and do not reproduce any board's questions.
These notes are free to read here, with no account.
What the Laisha app adds:
Free for ten days, then a subscription. We would rather write that here than have you find it later.