16  Consolidation and review

After working through this chapter you will be able to

• state in one paragraph what each chapter of this book established and how it depends on the others
• answer recall, application and interpretation questions across the whole syllabus and check your reasoning against worked answers
• solve problems that require combining methods from several chapters
• locate your own competence on a four-level rubric and name the next thing you cannot yet do
• choose a route through further study suited to a wet-lab, computational or clinical starting point

Fifteen chapters of techniques do not add up to a discipline until you can see the argument running through them. That argument is simple to state and takes a book to earn: biological questions become computable when biological objects are represented as data, and every subsequent decision — which alignment, which genetic code, which normalisation, which threshold — is a judgement about how much the representation has lost. A scoring matrix is a model of evolution. A count is not an abundance. A variant of uncertain significance is a statement about evidence, not about a variant. A predicted structure is a hypothesis with a confidence attached. The same idea returns in different vocabulary in every chapter.

A second thread runs alongside it, and it is the one most likely to save you: almost every expensive error in this book is silent. Prodigal exits zero and reports 995 genes where there are 524. A count matrix divided by a vector of library sizes recycles down the columns and returns plausible nonsense. Feature selection carried out before cross-validation reports perfect accuracy on randomly permuted labels. None of these raises an exception. What catches them is an independent check — a distribution you can predict, a reference you can compare against, a permutation that should destroy the signal.

This chapter consolidates both threads. The first section gives a tight synthesis of each of Chapters 1 to 15, written to be read in sequence so that the dependencies between them are visible. The second is a question bank of 66 graded questions split between recall, application and interpretation, with worked answers in collapsed callouts so that you can test yourself honestly before looking; wherever a chapter computed a number, the question is anchored to it. The chapter then supplies three integrative problems that cannot be solved from a single chapter, a self-assessment rubric with four competence levels, and three annotated learning pathways for readers arriving from a wet-lab, a programming or a clinical background.

Use it as a revision instrument rather than a replacement for the chapters. The syntheses are deliberately too compressed to learn from; they are written to tell you whether you have already learned the material, which is a different job. Where a synthesis contains a term you cannot define without looking it up, the honest conclusion is that the chapter it came from needs rereading, and the glossary in the appendices is the fastest route back.

import sys; sys.path.insert(0, "..")
from pybfx import book_data, style
style()

16.1 What each chapter established

Chapter 1: What bioinformatics is, and how to set up to do it

The opening chapter defined the discipline by what it produces: defensible biological statements extracted from large, noisy, machine-generated measurements, as distinct from computational biology, which produces methods and models. It argued that the scale is what forces computation, and did the arithmetic rather than asserting it — a single human genome at 30-fold coverage with 150 bp reads is about 620 million reads and 206 GB of uncompressed FASTQ, some 2,400 times the row limit of a spreadsheet. It then placed the field’s methodological eras against the data problem each one solved, and reduced almost every analysis to a five-stage shape in which each stage contains a decision that must be recorded in a script if the conclusion is to survive questioning months later. The chapter ended by building the environment the rest of the book uses and running a first real analysis: the SARS-CoV-2 reference genome at 29,903 bp and 38% GC, with long open reading frames concentrated on one strand. Everything that follows is a specialisation of that pattern — a question, a decision, a record.

Chapter 2: Molecular biology for people who analyse data

This chapter supplied the biology the algorithms encode, and taught it through the errors it prevents. Strandedness comes first, because a FASTA record stores one strand and says nothing about which, so the reverse complement — not the complement — is what every format expects. Translation never fails, which is why a plausible-looking protein is no evidence of a correct reading frame; the diagnostics are an initial methionine, a single terminal stop and no internal stops. Testing a naive open reading frame scan against the SARS-CoV-2 annotation made the point quantitatively: a 100-codon threshold recovers the long genes, misses four real ones and produces a chance hit on the empty strand, while lowering it to 38 codons finds all eleven genes at the cost of growing the output from 9 predictions to 76. The TP53 worked example — roughly 19 kb of genomic span, a 2,512 nt mRNA, an 1,182 nt coding sequence and a 393-residue protein identical to UniProt P04637 — establishes both the central dogma and the cheapest available test of a nucleotide-to-protein pipeline. Coordinates close the chapter, because 1-based inclusive and 0-based half-open systems never announce themselves.

Chapter 3: Biological databases and how to get data out of them

Here the book moved from biology to the records that represent it, and to the question of how much a record can be trusted. Resources divide into primary archives, which store what was submitted; curated knowledgebases, which store conclusions with evidence; and derived resources, which compute — and the kind of resource determines the sentence you are entitled to write about it. Accession discipline follows: INSDC mirroring makes GenBank, ENA and DDBJ one archive under one set of accessions; prefixes encode provenance, with NM_/NP_ curated and XM_/XP_ machine-predicted, and GCA_ an author’s assembly against GCF_ the RefSeq copy; and the digits after the full stop are a version, without which an accession does not identify data. Dissecting the UniProt entry for p53 showed what curation buys — 238 papers read, 1,120 cross-references and more than 1,500 features on a protein of 393 residues — and established the chapter’s most portable habit: read the evidence code before the annotation, because IDA and IMP are experiments while IEA and ISS are hypotheses. Programmatic access and recorded provenance — accession, version, date, checksum — make the rest of the book repeatable.

Chapter 4: Sequence alignment and what it tells you

Alignment is the field’s central algorithmic idea, and the chapter built it from the recursion up. Similarity is measured and homology is inferred, so no sequence is ever “40% homologous” to another; orthologues arise by speciation and usually keep their function, paralogues by duplication and need not. Dynamic programming gives an exact optimum under a stated scoring scheme, and the worked score matrix and traceback show what that optimum is made of. The chapter then demonstrated, rather than asserted, that the scoring scheme matters: the same distant globin pair aligns over 142 columns at 27.5% identity under BLOSUM62 or PAM250, but collapses to 101 columns at an apparently better 48.5% identity under PAM30, so a matrix chosen for the wrong evolutionary distance silently shortens an alignment instead of failing. BLAST trades optimality for speed, and its statistics carry the same lesson — the bit score of an alignment is 44.3 whatever database it was found in, while the E-value scales with database size. Multiple alignment of eight p53 orthologues localised conservation to the DNA-binding domain, and a bootstrapped tree showed which clades the data actually resolve and which, at 18% and 34% support, they do not.

Chapter 5: The command line as a working environment

The toolchain chapters begin with the observation that bioinformatics data are mostly line-oriented text files too large to open. The shell is presented as a programming environment organised around one idea — programs read a stream, write a stream, and pipes compose them — so that memory use stays constant however large the file. The text-processing core of grep, cut, sort, uniq, sed, awk and join is taught against the formats the book uses, and each idiom arrives with the specific silent error it prevents: sorting before uniq -c, anchoring a regular expression, checking row counts across a join, which is how the chapter discovers that 30 genes produce 31 joined rows because CDKN2A has two reviewed UniProt entries. An awk example that computes a perfectly plausible and entirely wrong GC content of 46.42%, through nothing worse than putting two statements in the wrong order, makes the chapter’s point about silent failure. Compression is treated as streaming rather than storage, with BGZF and tabix turning a compressed VCF into a randomly accessible one, and the chapter closes on exit status, set -euo pipefail, defensive scripts and SLURM job arrays.

Chapter 6: Python for bioinformatics

This chapter turns the reader from a user of tools into an author of them. It covers the language core — types, collections, control flow, functions, files — with every example drawn from sequence data, and insists on two habits: choose the collection from the question (a list for order, a dictionary for lookup, a set for membership) and convert coordinates in one tested place rather than in your head, because Python is 0-based and half-open while GFF3, VCF and GenBank are 1-based and closed. Biopython’s SeqRecord and feature.extract then remove the hand-arithmetic entirely, which the chapter verifies by translating the SARS-CoV-2 spike feature and recovering a 1,273-residue protein beginning MFVFLVLLPLVS. NumPy earns its place on a 20,000 × 150 quality matrix, where the mean Phred score falls from Q37.2 at the first cycle to Q23.7 at the last and crosses Q30 at cycle 111 — the kind of summary that decides a trimming parameter. pandas turns a GFF3 or a VCF into a table, and the chapter ends where professional code begins: two assertions are a test suite, and a module with a __main__ guard is a tool.

Chapter 7: R and Bioconductor

R enters where statistics and figures become the task. The chapter’s organising claim is that R’s design removes loops and index errors by applying operations to whole vectors — and that the same design produces the book’s most dangerous silent error, because recycling repeats a short operand down the columns of a matrix, so dividing a count matrix by a vector of library sizes returns plausible nonsense rather than an error. Factors, coercion and the tibble–data frame distinction are treated in the same spirit: each is a category of wrong answer rather than a language feature. The grammar of graphics is taught by debugging a figure, taking an unreadable linear-scale boxplot of counts spanning five orders of magnitude to a log scale and then to facets that reveal library protocol, not treatment, as the dominant source of variation in the pasilla data. Bioconductor’s value is located in its shared classes: Biostrings, GRanges with rtracklayer, and above all SummarizedExperiment, which makes the commonest serious bug in expression analysis — sample annotation detached from its counts — impossible to commit. A three-line differential test closes it: 1,547 genes pass an uncorrected threshold against roughly 407 expected by chance, and 195 survive Benjamini–Hochberg correction.

Chapter 8: From assembled sequence to annotated genome

With sequence, tools and languages in hand, the book turns to interpretation. An assembly is a hypothesis, and the chapter checks contiguity, completeness and contamination before annotating anything, noting that N50 measures contiguity rather than correctness and rewards over-aggressive assembly. Structural annotation of Mycoplasmoides genitalium then delivers the chapter’s central lesson about parameters: under the bacterial default, genetic code 11, Prodigal predicts 995 genes where the correct code 4 predicts 527, and nothing in the output signals the error. Comparing the prediction with the RefSeq annotation of the same assembly shows what ab initio prediction is and is not good at — 464 genes agree on both coordinates, 37 share a stop but not a start, 23 are found only by RefSeq and 26 only by Prodigal — and explains the asymmetry: a stop codon is a consequence of the reading frame while a start is a prediction. Eukaryotic prediction needs external evidence because introns make the problem underdetermined. Functional annotation is inference by similarity with an invisible inference chain, which is why the chapter ends by reading evidence codes: all 985 Gene Ontology annotations in this genome are IEA, and about a tenth of the proteome remains unannotated after thirty years.

Chapter 9: Next-generation sequencing data from reads to results

This chapter connects an instrument to a result and insists that each stage is a hypothesis rather than a view of the same thing: a read is an observation, an alignment an inference about origin, a variant call an inference about genotype. It relates each chemistry to the errors it makes, separates library artefacts from biology, and builds a quality report from the bundled reads rather than trusting a tool’s traffic lights. Trimming is quantified rather than recommended: at Q25 the data loses 11% of its bases for no gain in an alignment-based analysis. Alignment introduces SAM in detail — FLAG, CIGAR and MAPQ — and the sharpest demonstration uses simulated reads whose true origin is known: 25 of 40,000 reads are placed in the wrong location, 38 carry MAPQ 0 in the MgPa adhesin repeats, and MAPQ 0 means “equally good elsewhere” rather than “bad read”. The RNA-seq half runs a full DESeq2 analysis of the pasilla experiment, recovering 1,108 genes at an adjusted p-value below 0.05 and pasilla itself down about 3.7-fold, and ends on what a 5% false discovery rate is a property of: the list, not the gene.

Chapter 10: Protein structure and what you can compute from it

Structure is where sequence becomes mechanism, and the chapter treats a structure as a measurement with error bars. Judging one means reading resolution, R factors, B factors and missing residues — 1TUP at 2.20 Å reports no R-free because it was deposited in 1995, and 72 of its residues are modelled nowhere, recorded in REMARK 465 where code that reads only ATOM lines will never see them. Bio.PDB’s structure–model–chain–residue–atom hierarchy makes geometry a few lines of code, and the geometry is biology: the p53 zinc is held tetrahedrally by Cys176, His179, Cys238 and Cys242 at 1.99–2.33 Å without touching DNA, while only chain B is bound in the major groove, through twelve residues within 4 Å of it. Superposition supplies the chapter’s most transferable number — three copies of the same domain in one crystal differ by 0.78–1.03 Å, which is the scale against which any other RMSD must be read. The AlphaFold model reproduces that domain to 0.73 Å over 194 Cα atoms and scores p53’s disordered regions below 50, which is a correct prediction rather than a failure; high predicted aligned error between two confident domains means their relative arrangement is not predicted at all.

Chapter 11: Biological networks and pathway analysis

Having treated genes one at a time, the book turns to their relations, and spends most of its effort on what an edge actually asserts. In the bundled 30-gene STRING slice every edge has text-mining support, 176 of 210 have experimental support and only 85 appear in a curated database, so a functional association is not a physical interaction and filtering by evidence channel is usually the right first move. The chapter then measures the network — 30 nodes, 210 edges, density 0.483, clustering coefficient 0.646 — and immediately disowns the measurements as properties of a query-driven selection rather than of the interactome. The same scepticism applies to hubs: TP53 has degree 29 and five times the betweenness of any other gene, which in a network built by querying for p53-associated genes is no discovery, and degree correlates with annotation count at Spearman 0.746, so hub status partly measures study intensity. Community detection always returns communities, so low modularity and disagreement between algorithms are the signal. Over-representation analysis is then built from one hypergeometric test per term, corrected with Benjamini–Hochberg, and the closing demonstration is that matching the background universe for study intensity cuts the significant-term count from 511 to 121.

Chapter 12: Statistics and machine learning on omics data

This chapter supplies the statistical reasoning the rest of the book consumes, and does it by breaking things deliberately. It distinguishes biological from technical replication, establishes that an omics matrix has far more features than samples — 14,599 genes against seven pasilla samples — and shows what that guarantees: some feature always separates any two groups. Normalisation is demonstrated rather than described, with the median log ratio between two control samples falling from 1.443, the log ratio of their library sizes, to 0.056 after conversion to counts per million. Batch effects are shown imitating biology: principal components place treatment on PC1 at 39.7% of variance and library protocol on PC2 at 31.8%, and in a deliberately confounded contrast 22.3% of the significant genes are also significant in a pure protocol comparison against 7.8% expected. Multiple testing is calibrated against simulation — 10,000 null tests give 501 results at p < 0.05 and none after correction — and the most important experiment permutes the labels: selecting features before cross-validation reports perfect accuracy in 100 of 100 permutations of meaningless labels, while selecting inside each fold reports 0.384, below the majority-class baseline, which is the truthful answer.

Chapter 13: Bioinformatics in disease and healthcare

Here the book’s methods meet the clinic, and the chapter’s discipline is to separate claims by the evidence behind them. Germline and somatic variants differ in calling strategy, expected allele fraction and consent obligations; drivers are separated from passengers by comparison against a background mutation rate that varies by orders of magnitude; mutational signatures describe the chemistry that produced a mutation and not its consequence. The bundled ClinVar slice of 3,888 TP53 records supplies the clinical reality: about a quarter are variants of uncertain significance, only 191 carry expert-panel review and 192 carry no assertion criteria at all, so review status matters more than the label. Two computed results carry the chapter. Truncating variants are classified pathogenic 81.4% of the time against 16.9% for missense variants, which is the empirical content of the ACMG/AMP criterion PVS1. And 94.1% of pathogenic missense variants fall in the DNA-binding domain at residues 102–292, against 26.1% of benign ones and the 48.6% expected from its length. The chapter closes on the distinction between pathogenic and actionable, on pharmacogenomic guidelines, and on pathogen surveillance, which is limited by sampling and reporting far more often than by sequencing.

Chapter 14: Working responsibly with genomic data

This chapter argues that the obligations attached to human genomic data follow from four properties of the data itself: it identifies, it cannot be revoked, it implicates relatives who never consented, and it predicts weakly but is read strongly. Two simulations make the first property quantitative. Twenty-five independent common markers make every one of 10,000 simulated individuals unique, and the random match probability falls below one in eight billion at 30 markers, which is why a variant file is identifiable data whatever has been stripped from its header. A reimplementation of the allele-frequency membership statistic then detects 90% of 200 pool contributors from 5,000 markers with no false positives, showing that published summary statistics are not automatically safe and that the attack strengthens as panels densify. From that mechanical basis the chapter works through consent models and the problem of enforcing them in code, controlled access and the five safes, the structure of GDPR, HIPAA, Nigeria’s data protection regime and the Nagoya Protocol, the clinical cost of the European-ancestry bias in reference data, and the design decisions — benefit sharing, authorship, governance — that separate collaboration from extraction. Its practical residue is a checklist applied before a project starts.

Chapter 15: A complete project, start to finish

The capstone assembles everything into one analysis: how much of a curated bacterial annotation does ab initio prediction recover? It builds the infrastructure first — read-only raw data, disposable derived directories, a pinned environment, a Git history whose messages record decisions, and a Makefile that rebuilds only what is stale — and then runs an analysis that fails. Prodigal’s default translation table treats TGA as a stop codon, which 71.6% of M. genitalium genes contain internally, so the first run returns 995 predictions with a median length of 336 bp against 856 bp in the reference, recovering 96.6% of TGA-free genes but only 52.5% of the rest. Diagnosing that against an independent annotation, rather than from an error message, is the chapter’s central lesson. Under the correct code the comparison recovers 95.6% of reference coding sequences at 95.1% precision with 92.6% of starts identical, and inspecting the residual disagreements shows them concentrated in very short genes and in the MgPa adhesin repeat family, where neither method is simply right. It closes on write-up: the metric belongs in the abstract, the limitations are part of the result, and a reviewer’s checklist is cheaper applied to yourself.

16.2 Question bank: foundations and data

Each question is labelled by the kind of thinking it requires. Recall questions test whether a definition or fact is available to you. Application questions require you to use a method on a specific case. Interpretation questions give you a result and ask what it means or whether it is trustworthy; these are the ones that matter most in practice and the ones most often skipped. Attempt each before opening the answers.

Chapter 1: What bioinformatics is, and how to set up to do it

1. (recall) State the distinction Chapter 1 draws between bioinformatics and computational biology, in terms of what each produces.
2. (recall) Name the five stages common to almost every sequence analysis, and state what all five have in common.
3. (application) A 30-fold human genome sequenced with 150 bp reads produces about 620 million reads and 206 GB of uncompressed FASTQ. Show the arithmetic, and explain why the comparison with a spreadsheet’s row limit is a statement about operations as much as about size.
4. (interpretation) A published analysis reports a result but gives neither software versions nor parameters. State what specifically cannot be checked, and why “the method is standard” is not an adequate defence.

Chapter 2: Molecular biology for people who analyse data

5. (recall) Define reading frame, and state how many exist for a double-stranded DNA sequence.
6. (recall) Why is a plausible-looking translated protein not evidence that the reading frame was correct, and what three features should you check instead?
7. (application) A single-base deletion occurs early in the TP53 coding sequence. Predict the consequence for the protein product and justify it from the triplet nature of translation.
8. (interpretation) Scanning the SARS-CoV-2 genome for open reading frames of at least 100 codons recovers the long genes but misses four annotated ones and reports a hit on the strand that carries no genes. Lowering the threshold to 38 codons finds all eleven genes but raises the output from 9 predictions to 76. State what this trade-off is called and how you would choose a threshold.

Chapter 3: Biological databases and how to get data out of them

9. (recall) Distinguish a primary archive, a curated knowledgebase and a derived resource, naming one of each.
10. (recall) What do the digits after the full stop in NM_000546.6 specify, and what is versioned differently in an Ensembl identifier such as ENST00000269305.9?
11. (application) You need the canonical human p53 sequence and the boundaries of its DNA-binding domain. State which resource you would query, where in the record the boundaries sit, and what you would record as provenance.
12. (interpretation) Two functional annotations on the same protein carry the evidence codes IDA and IEA. State what each means, which you would build a conclusion on, and what you would check before relying on the weaker one.

Chapter 4: Sequence alignment and what it tells you

13. (recall) State the difference between global and local alignment, name an algorithm for each, and give the biological question that selects each one.
14. (recall) Explain what an affine gap penalty models that a linear penalty does not.
15. (application) The same distant globin pair aligns over 142 columns at 27.5% identity under BLOSUM62 and over 101 columns at 48.5% identity under PAM30. Explain why the higher identity is the worse result, and state which matrix you would report.
16. (interpretation) A BLAST search of myoglobin against a small protein database returns human haemoglobin alpha at 25.7% identity over 148 residues with an E-value of 9.18 × 10⁻¹⁰, and several p53 entries matching over 8 to 27 residues with E-values between 0.47 and 9.7. State which hits are evidence of homology, and why percentage identity is the weaker guide.

Answers: Chapters 1 to 4

1. Bioinformatics converts large, noisy, machine-generated measurements into defensible biological statements about a particular sample; computational biology produces the methods and models that make such conversions possible. The distinction is in the product, not in the tools, and it determines how the work is judged — correctness, reproducibility and relevance to the sample in hand for the first, generality and correctness of the method for the second.

2. Acquiring or retrieving the data; assessing its quality; reducing or transforming it; comparing it with a reference or between groups; and interpreting the comparison. What they have in common is that every stage contains a decision — a threshold, a reference, a parameter — which must be recorded in a script, because the question asked months later is whether the conclusion survives changing it.

3. A 3.1 Gb genome at 30-fold coverage is about 93 Gb of sequence; divided by 150 bp that is roughly 620 million reads, and at about 332 bytes per FASTQ record that is some 206 GB. The spreadsheet comparison matters less for the row count than for the operations: aligning each read to a 3.1 Gb reference, modelling base distributions under an error profile, and testing 20,000 genes with false discovery control are not row-by-row arithmetic, so no amount of spreadsheet capacity would help.

4. Without versions you cannot know whether the tool’s defaults, algorithms or bug fixes differ from yours, and without parameters you cannot reproduce the numbers at all; the result is unverifiable rather than merely inconvenient. “Standard” is no defence because defaults change between versions and because the parameter that mattered most in this book — the translation table in Chapters 8 and 15 — was itself a default.

5. A reading frame is one of the ways a nucleotide sequence can be partitioned into consecutive triplets, fixed by the position of the first base read. There are six for double-stranded DNA: three on each strand.

6. Because translation never fails — any sequence of bases yields a sequence of amino acids, so the output always looks like a protein. Check for an initial methionine, a single stop at the end, and no internal stops; a frame-shifted translation typically fails the last two.

7. The deletion shifts the frame for every downstream codon, so the sequence after it is wrong and a premature stop appears quickly. In the chapter’s worked example ten correct residues are followed by a short stretch belonging to no p53 and then a stop — 13 amino acids in place of 393, from one base in a 1,182 nt coding sequence. The transcript may also be destroyed by nonsense-mediated decay.

8. It is the sensitivity–precision trade-off: lowering the threshold recovers the short genes (sensitivity) and admits many chance frames and alternative starts (precision), since roughly six of every seven predictions at 38 codons are not new genes. Choose the threshold from what the prediction is for — a high threshold when you will inspect every hit by hand, a low one when the predictions feed a downstream filter such as homology search — and state the threshold with the result.

9. A primary archive such as GenBank stores records as submitted, with the submitter responsible for their content; a curated knowledgebase such as UniProtKB/Swiss-Prot or RefSeq records conclusions with evidence attached; a derived resource computes something from the others, such as a genome browser’s alignment tracks or a structure prediction database.

10. In an INSDC or RefSeq nucleotide accession the version increments when the sequence changes, so NM_000546.5 and NM_000546.6 may be different strings of bases. Ensembl versions genes and transcripts separately from the assembly and increments the suffix when the model changes, so the same genomic sequence can carry a new transcript version.

11. Query UniProtKB for P04637. The domain boundaries are in the feature table of the flat file, as the DNA_BIND 102..292 line, with its own evidence attribution. Record the accession, the entry version or release date, the date of retrieval and a checksum of the downloaded file, because UniProt annotations change.

12. IDA is inferred from direct assay — an experiment on this protein — while IEA is inferred from electronic annotation, transferred automatically from another record without curator review. Build on the IDA annotation. Before relying on the IEA one, find what it was transferred from, check that the similarity covers the functionally relevant region, and look for direct experimental evidence; an electronic annotation is a hypothesis about function, not an observation of it.

13. Global alignment aligns two sequences along their entire length (Needleman–Wunsch) and suits sequences expected to correspond throughout, such as two orthologues of the same protein. Local alignment finds the best-scoring segment pair (Smith–Waterman) and suits the search for a shared domain in otherwise unrelated sequences.

14. That an insertion or deletion is a single evolutionary event of some length, rather than a series of independent single-residue events. A linear penalty charges the same for each gap position and therefore scatters many short gaps where one long gap belongs; an affine penalty charges a large opening cost and a small extension cost, which recovers the single event.

15. PAM30 is built for sequences that have barely diverged, so it penalises the substitutions that actually occur between distant homologues and the aligner responds by truncating the alignment to a short, well-matching core. The identity rises because the denominator shrank, not because the sequences are more similar. Report the BLOSUM62 or PAM250 result — 142 columns at 27.5% — and report identity with the aligned length attached, since identity alone is uninterpretable.

16. Only the haemoglobin alpha hit. The E-value accounts for alignment length and database size, so 9.18 × 10⁻¹⁰ over 148 residues is strong evidence of homology even at 25.7% identity, which sits in the twilight zone where real homologues and convincing accidents coexist. The p53 matches over 8 to 27 residues have E-values close to or above 1, meaning roughly one such hit is expected per search by chance. Percentage identity ignores both the length of the alignment and the size of the search space.

16.3 Question bank: the toolchain

Chapter 5: The command line as a working environment

17. (recall) Explain what a pipe does, and why it allows processing of files much larger than available memory.
18. (recall) State the difference between > and >>, and the specific data loss that confusing them causes.
19. (application) Write a one-line pipeline that counts the sequences in a FASTA file, and explain why counting lines beginning with > is correct whereas counting all lines is not.
20. (interpretation) An awk program reports the GC content of a read file as 46.42%, a plausible figure, when the true value is near 31.7%. The program removes the newline before accumulating base counts but after accumulating the total length. Explain the error, and state the general lesson about this class of bug.

Chapter 6: Python for bioinformatics

21. (recall) State what SeqIO.parse returns and how it differs from SeqIO.read.
22. (recall) Python slices are 0-based and half-open; GFF3, VCF and GenBank coordinates are 1-based and inclusive. Give the Python slice for the GFF3 feature 686..1828, and state the rule.
23. (application) Write a function returning the GC content of a sequence that handles an empty sequence without raising. State what it should return and justify the choice.
24. (interpretation) The mean per-cycle quality of the bundled reads falls from Q37.2 to Q23.7 and crosses Q30 at cycle 111, while the tenth percentile crosses it at cycle 83. State which of these two numbers should inform a trimming decision and why.

Chapter 7: R and Bioconductor

25. (recall) State what recycling is in R, and the direction in which it operates on a matrix.
26. (recall) What three behaviours does a tibble change relative to a data frame, and what class of error does each prevent?
27. (application) You have a counts matrix and a separate sample table. Describe the checks to perform before any analysis, and the single assertion that must hold between them.
28. (application) Explain how pivot_longer changes the shape of a 14,599-gene by 7-sample count matrix, why ggplot2 needs that shape, and why Bioconductor assays need the opposite one.
29. (interpretation) A boxplot of per-sample gene counts on a linear axis shows flattened boxes and a column of outliers. State the diagnosis and the single change that fixes it, and say why changing the geometry would not.
30. (interpretation) Faceting the pasilla count distributions by library protocol shows medians differing by roughly a factor of two between panels, while the four distributions within each panel are closely matched. State what this identifies, and what it implies for any differential expression model fitted to these samples.

Answers: Chapters 5 to 7

17. A pipe connects one program’s standard output to the next program’s standard input, so records stream through the chain one at a time. Because no stage holds the whole file, memory use stays roughly constant regardless of file size.

18. > truncates the target file before writing; >> appends. Confusing them destroys the previous contents irrecoverably, which is why raw data directories are made read-only with chmod.

19. grep -c "^>" sequences.fasta. Only header lines start with >, and each record has exactly one header, whereas sequence lines are wrapped at an arbitrary width, so a line count measures formatting rather than content. The ^ anchor matters: without it the count would include any line containing >.

20. The length accumulator sees the line with its trailing newline and the base counts do not, so every line contributes one phantom base to the denominator; with 150 bp lines the inflation is about 0.7% of the total and the GC percentage shifts by a plausible amount. The lesson is that this class of error produces a number in the right range rather than a crash, so a pipeline needs an external check — here, the genome’s own known GC content of 31.7% — rather than a plausibility judgement.

21. SeqIO.parse returns an iterator over all records in a file and suits a file with any number of sequences; SeqIO.read returns the single record in a file and raises an error if there is not exactly one. Use read when “exactly one record” is a requirement you want enforced.

22. genome[685:1828]. Subtract one from the start and leave the end unchanged; the half-open upper bound and the inclusive end coincide numerically, which is why only the start appears to move.

23. It should return a documented sentinel such as None or float("nan"), not zero and not an exception, because an empty sequence has no GC content and 0% would enter a downstream average as though it were a measurement.

def gc_fraction_safe(sequence):
    """GC content as a fraction, or None for an empty sequence."""
    if len(sequence) == 0:
        return None
    gc = sum(1 for base in sequence.upper() if base in "GC")
    return gc / len(sequence)

24. The tenth percentile, at cycle 83. The mean describes the typical read, but mismatches are contributed by the worst reads, and a cycle at which a tenth of reads have already fallen below the accuracy threshold is adding more error than information. Trimming decisions are made from the lower edge of the distribution, not from its centre.

25. Recycling repeats the shorter operand until it matches the length of the longer one, silently and without warning when the lengths divide evenly. A matrix in R is a vector with a dimension attribute filled column by column, so a short vector recycles down the columns — which is why dividing a gene-by-sample count matrix by a vector of per-sample library sizes scales genes rather than samples and returns a plausible, wrong answer.

26. No partial matching of column names, which prevents a typo silently selecting the wrong column; no dropping to a vector when a single column is selected, which prevents code that expects a table receiving something else; and no silent conversion of character columns to factors, which prevents the whole class of factor-related errors.

27. Check that the number of samples equals the number of matrix columns, that the two identifier sets are identical, and that they are in the same order; then reorder the sample table by the matrix column names rather than assuming correspondence. The assertion that must hold is that sample-table row names equal matrix column names element by element — or, better, build a SummarizedExperiment, which enforces it permanently.

28. pivot_longer turns the 14,599 × 7 wide matrix into roughly 102,000 rows of three columns — gene, sample, count — one row per observation. ggplot2 needs this because every aesthetic maps to a column, so “sample” must be a variable rather than a set of column names. Bioconductor assays need the wide form because the matrix is the computational object: linear algebra, normalisation and model fitting all operate on genes × samples.

29. The counts span five orders of magnitude, so a linear axis gives almost all its height to a few thousand very large values and compresses everything else onto the axis. Change the scale to logarithmic; the geometry is not the problem, and replacing the boxplot with a violin or a histogram would produce an equally unreadable figure on the same axis.

30. It identifies library protocol, not treatment, as the dominant source of variation — a batch effect of larger magnitude than the biological effect. Since protocol is recorded and is not perfectly confounded with condition in this dataset, it must appear in the design formula as a covariate, so that the condition term is estimated after accounting for it; fitting condition alone would attribute protocol variation to treatment and inflate the dispersion estimates.

16.4 Question bank: genomes, reads and structures

Chapter 8: From assembled sequence to annotated genome

31. (recall) Distinguish structural from functional annotation, and name one tool used for each.
32. (recall) State what the nine columns of a GFF3 record contain, and which attributes make a feature’s relationship to another explicit.
33. (application) You have a new bacterial assembly and no annotation. Give the order of steps to a first-pass annotation, and say at which step the translation table must be settled and why.
34. (interpretation) Comparing Prodigal with the RefSeq annotation of the same genome gives 464 genes identical at both ends, 37 sharing a stop but not a start, 23 found only by RefSeq and 26 only by Prodigal. Explain why the disagreements concentrate in start codons rather than in gene presence, and what the 23 and the 26 mostly turn out to be.

Chapter 9: Next-generation sequencing data from reads to results

35. (recall) State what a Phred score of 30 means, and how the score is encoded in a FASTQ quality line.
36. (recall) Name the eleven mandatory SAM fields’ three most frequently misread members — FLAG, CIGAR and MAPQ — and state what each encodes.
37. (application) Given paired FASTQ files, give the commands that produce a sorted, indexed BAM, and state why coordinate sorting must precede indexing.
38. (interpretation) samtools flagstat on the bundled simulated alignment reports 40,000 reads, 100% mapped and 39,968 properly paired. Say what each figure would mean for a real bacterial isolate, and which of the three you would be most suspicious of here.
39. (interpretation) Thirty-eight of 40,000 reads carry MAPQ 0, clustered around two positions, and 25 reads are placed somewhere other than where they were simulated from. Explain what MAPQ 0 means, why the two observations are related, and what filtering on MAPQ 20 would do to a variant call in that region.
40. (interpretation) The DESeq2 analysis of the pasilla data returns 1,108 genes at an adjusted p-value below 0.05, and recovers pasilla itself as down about 3.7-fold. State precisely what the 5% false discovery rate promises about this list, and what it promises about any one gene in it.

Chapter 10: Protein structure and what you can compute from it

41. (recall) State what the B-factor column reports in a crystal structure, and what occupies the same column in an AlphaFold model.
42. (recall) Where in a PDB file are residues that were not modelled recorded, and why will code that reads only ATOM records miss them?
43. (application) Describe how you would determine which residues of a protein chain contact a bound DNA molecule, including the atom selection and the cut-off, and what you would exclude from both sets.
44. (interpretation) Three copies of the p53 DNA-binding domain in one crystal differ by 0.78–1.03 Å over about 195 Cα atoms, and the AlphaFold model matches one of them to 0.73 Å over 194 Cα. State what the first range is for, and what the second number does and does not establish.
45. (interpretation) An AlphaFold model has two regions with pLDDT above 90 and a high predicted aligned error between them. State what each measure reports and which statements the model does and does not support.
46. (interpretation) The radius of gyration of the p53 DNA-binding domain alone is 16.6 Å, while that of the full-length AlphaFold model is 33.5 Å. Explain what the difference reflects, and why the full-length figure should not be quoted as a property of the protein.

Answers: Chapters 8 to 10

31. Structural annotation locates features — genes, exons, RNA genes, repeats — in an assembly; Prodigal performs it for bacteria and AUGUSTUS for eukaryotes. Functional annotation assigns biological meaning to those features; domain search against Pfam with HMMER, or InterProScan, performs it.

32. Sequence identifier, source, feature type, start, end, score, strand, phase and attributes. The ID and Parent attributes in the ninth column express the hierarchy, making a CDS a child of a transcript and a transcript a child of a gene; a spliced feature occupies several lines sharing one ID.

33. Assess the assembly — contiguity, completeness against a conserved single-copy gene set, contamination — then mask repeats, then predict coding sequences, then predict RNA genes with specialised tools, then transfer functional annotation by domain and orthology search, and finally assess the annotation before release. The translation table must be settled before gene prediction, because it determines which codons terminate an open reading frame: on M. genitalium the bacterial default produces 995 predicted genes against 527 under the correct code 4, and nothing in the output signals the error.

34. Because a stop codon is a consequence of the reading frame — once the frame and the 3’ end are right, both methods must agree — whereas a start codon is a prediction made from a ribosome binding site model among several in-frame candidates within the first 60 bp. Both unmatched groups turn out to be the same handful of loci seen from two sides. Nineteen of the 23 RefSeq-only genes coincide with annotated pseudogenes, which homology evidence supports but which carry a frameshift or internal stop that an ab initio model will not call; and 20 of those same 23 do overlap a Prodigal prediction, so they are boundary disagreements rather than genes the predictor failed to see at all. Correspondingly, 19 of the 26 Prodigal-only predictions overlap one of the 20 annotated pseudogene loci. Only three reference genes are absent from the prediction entirely, and the shortest of them, a 147 bp ribosomal protein gene, shows where the length sensitivity of ab initio prediction runs out.

35. A Phred score of 30 means an estimated one-in-1,000 probability that the base call is wrong, since \(Q = -10\log_{10}P\). It is stored as a single ASCII character whose code point is \(Q + 33\), so \(Q = 0\) is ! and \(Q = 40\) is I; files from Illumina pipelines before version 1.8 used an offset of 64 instead, which makes the same character mean a score 31 points higher.

36. FLAG is a bit field recording pairing, strand, mate strand, mate mapping, duplicate status and whether the alignment is secondary or supplementary. CIGAR describes how the read aligns — matches, insertions, deletions, soft clipping — and therefore determines where the alignment ends, which the read length alone does not. MAPQ is the Phred-scaled confidence that the read was placed in the right location.

37. Align with bwa mem piped into samtools view to produce BAM, sort with samtools sort, then index with samtools index. The index maps genomic coordinates to file offsets so a region query can seek directly to the right block; that mapping only exists if the records are stored in coordinate order.

38. For a real bacterial isolate, 95–99% mapped is expected and anything below about 70% calls for taxonomic classification of the unmapped reads; the properly paired figure adds the constraint that mates are in the expected orientation and separation, and a low value points at adapter or insert-size problems. All three figures here are suspicious because they are too good: 100% mapped and 99.92% properly paired are the signature of reads simulated from the same reference, as is the insert size of exactly 349 bp with zero standard deviation.

39. MAPQ 0 means the read aligns equally well in at least two places and the aligner reported one of them arbitrarily — it is a statement about the reference, not about read quality. The two observations are the same phenomenon: the ambiguous reads sit in copies of the MgPa adhesin repeat family, and a read whose true origin was one copy can be placed in the other. Filtering at MAPQ 20 discards those alignments, which is usually correct because a variant supported only by ambiguously placed reads is an artefact; the cost is that genuine variants in the repeat become invisible rather than wrong, and that difference must be stated in the methods.

40. It promises that, among the 1,108 genes on the list, about 5% are expected to be false discoveries. It promises nothing about any individual gene — a gene with an adjusted p-value just under the threshold is not “95% likely to be real” — and nothing at all about effect size or biological importance, which is why the fold change and the mean expression must be read alongside it. Recovering the knocked-down gene itself, without having been told to look for it, is a positive control on the whole pipeline rather than part of the statistical claim.

41. In a crystal structure the B-factor reports positional uncertainty or thermal motion for each atom. In an AlphaFold model the same column holds pLDDT, the per-residue confidence on a 0–100 scale, which is not a temperature factor and must not be interpreted as one.

42. In REMARK 465, not in the coordinate records. Code that parses only ATOM lines sees a continuous chain where there is a gap — 72 residues of 1TUP are modelled nowhere, all at chain termini — so distances computed across such a gap are meaningless and a sequence read from the coordinates is incomplete.

43. Take all atoms of standard residues in the protein chain and all atoms of standard residues in the DNA chains, compute the minimum inter-atomic distance from each protein residue to the DNA, and keep residues below a cut-off of about 4 Å, which is the range of hydrogen bonds and close van der Waals contacts. Exclude waters, ions and other heteroatoms from both sets: leaving ordered water molecules in the DNA set inflates the contact list with residues that touch solvent rather than nucleic acid.

44. The 0.78–1.03 Å range is a calibration: it is how much the same domain differs from itself within one crystal, so it is the noise floor against which any other RMSD must be read. The AlphaFold figure of 0.73 Å over 194 Cα establishes that the prediction of that domain is as close to the experiment as the experimental copies are to each other, which is as good as the comparison can resolve. It does not establish that the rest of the model is accurate, nor anything about the domain’s orientation relative to the other parts of the protein.

45. pLDDT measures local confidence — whether the backbone around each residue is placed correctly — and predicted aligned error measures confidence in the relative position of one part of the model with respect to another. Two confident domains with high error between them support statements about each domain’s internal geometry, such as the shape of a binding pocket, and support no statement about their arrangement, so inter-domain distances, overall shape and any docking that spans the two are unsupported.

46. The difference reflects the disordered regions tethering the two confident domains: 16.6 Å is what a compact 150–200 residue domain should give, while 33.5 Å describes a pair of domains drawn in one arbitrary arrangement on flexible linkers. Quoting it as a property of the protein implies a compactness the model does not assert — the arrangement is exactly what the predicted aligned error says is not predicted — and in solution the molecule samples an ensemble rather than this one conformation.

16.5 Question bank: systems, data science and application

Chapter 11: Biological networks and pathway analysis

47. (recall) State what an edge in a STRING network asserts, and how that differs from an edge in a curated physical interaction database.
48. (recall) Define degree centrality and betweenness centrality, and state what each identifies.
49. (application) Given 80 genes from a differential expression analysis, describe how you would test for enrichment of a biological process, and state what the background universe must be.
50. (interpretation) In the bundled network TP53 has degree 29 of a possible 29 and five times the betweenness of the next gene, and degree correlates with Gene Ontology annotation count at Spearman 0.746. State the two separate reasons this hub status is not a discovery.

Chapter 12: Statistics and machine learning on omics data

51. (recall) State what the false discovery rate controls, and how it differs from the family-wise error rate.
52. (application) You will train a classifier on seven samples and 14,599 features. Describe a validation scheme that does not leak information, name the step most often performed incorrectly, and give the check that detects the leak.
53. (interpretation) Before normalisation the median log ratio between two untreated pasilla samples is 1.443; after conversion to counts per million it is 0.056. Explain what the first number is, and what would have happened to a differential expression analysis run on the raw counts.
54. (interpretation) A classifier on a test set with 2.4% positives has an ROC area of 0.858 but an average precision of 0.407 against a no-skill baseline of 0.024, with precision 0.706 and recall 0.308. State which summary is honest for this problem and what the model is actually doing.

Chapter 13: Bioinformatics in disease and healthcare

55. (recall) State what a variant of uncertain significance is, and why its existence is a statement about evidence rather than about the variant.
56. (application) A patient’s exome contains a rare missense variant in a gene plausibly related to their phenotype. Give the steps you would take to assess it, in order, before any clinical statement is made.
57. (interpretation) In the bundled ClinVar slice, truncating TP53 variants are classified pathogenic 81.4% of the time against 16.9% of missense variants. Name the ACMG/AMP criterion this quantifies, and explain why the missense figure is not evidence that missense variants are usually benign.
58. (interpretation) 94.1% of pathogenic missense variants in TP53 fall in the DNA-binding domain, against 26.1% of benign ones and 48.6% expected from its length. State what this supports, which ACMG/AMP criterion it licenses, and what it does not tell you about any single variant.

Chapter 14: Working responsibly with genomic data

59. (recall) State why a genotype cannot be anonymised in the way a name can be removed, and give the quantitative basis for the claim.
60. (recall) Name the five dimensions of the five safes framework.
61. (application) You are asked to publish per-variant allele counts for cases and controls in a study of a stigmatised condition. State the attack this enables, the three inputs the attacker needs, and two modifications that preserve most of the scientific value.
62. (interpretation) A polygenic score validated in a European-ancestry cohort is proposed for clinical use in West Africa. State the mechanism by which its accuracy degrades, and what evidence would be required before deployment.

Chapter 15: A complete project, start to finish

63. (recall) State what belongs in a version-controlled repository for an analysis project and what must never be committed.
64. (recall) Explain what a Makefile records that a numbered list of scripts does not, and what make clean && make tests.
65. (application) Describe the test that establishes whether a project is reproducible, in the order the commands would be run.
66. (interpretation) Under genetic code 11 the capstone’s first run recovered 96.6% of the 149 reference genes with no internal TGA but only 52.5% of the 375 that contain one. Explain why this single comparison identifies the fault, and name the two external checks that made the failure visible when the tool itself reported none.

Answers: Chapters 11 to 15

47. A STRING edge asserts a functional association supported by some combination of evidence channels — co-expression, co-occurrence, text mining, databases, experiments — combined into a confidence score. A curated physical interaction edge asserts that two proteins were observed to bind. In the bundled subnetwork every edge has text-mining support, 176 of 210 have experimental support and only 85 appear in a curated database, so treating the network as a map of physical contacts would overstate two-thirds of it.

48. Degree centrality is the number of edges at a node and identifies hubs. Betweenness centrality is the fraction of shortest paths passing through a node and identifies bottlenecks connecting otherwise separate parts of the network. A node can have high betweenness and modest degree, which is often the more interesting case.

49. Map the genes to the identifier space of the annotation source, test each term for over-representation with a hypergeometric test, and correct across terms with Benjamini–Hochberg. The background universe must be the set of genes the experiment could have detected — the genes with non-zero counts in the assay, not the whole proteome — because testing against an inflated universe credits you for the enrichment of expressed genes. The chapter’s demonstration of this is blunt: matching the universe for study intensity cut the significant-term count from 511 to 121.

50. First, the network was built by querying a database for genes associated with p53 and the cell cycle, so TP53 is adjacent to every other node by construction; the topology reproduces the selection criterion. Second, degree in a literature-derived network partly measures how much a protein has been studied, which the Spearman correlation of 0.746 with annotation count makes explicit — these genes carry a median of 54.5 Gene Ontology terms against a genome-wide median of nine.

51. The false discovery rate is the expected proportion of false positives among the tests declared significant; the family-wise error rate is the probability of making even one false positive anywhere. FDR is the appropriate choice when thousands of hypotheses are tested and a known, small proportion of false discoveries is tolerable, which is why it leaves 511 enriched terms where Bonferroni leaves 82.

52. Hold out a test set before anything else, or use nested cross-validation, with normalisation, feature selection and hyperparameter tuning performed inside each training fold only. The step most often performed incorrectly is feature selection on the whole dataset before splitting. The check is label permutation: with labels randomly shuffled there is nothing to learn, so a leaky pipeline that reports a mean accuracy of 1.000 in 100 of 100 permutations has been caught, while an honest one reports a value scattered around or below the majority-class baseline — 0.384 against a baseline of 0.571 in the chapter’s demonstration.

53. It is the log ratio of the two libraries’ sizes showing up as an apparent expression difference: most genes do not change between two control samples, so the median log ratio measures sequencing depth rather than biology. Run on raw counts, a differential expression analysis would have called a large, systematic and entirely artefactual set of changes in the direction of the deeper library, and the effect would have been reproducible, which is what makes it dangerous.

54. The precision–recall summary is the honest one. With 2.4% positives the ROC area is flattered because the false positive rate is diluted by a large number of true negatives, so 0.858 sounds respectable while the model misses nearly 70% of the positives. What it is actually doing is being right when it commits — precision 0.706 — at the cost of committing rarely; whether that is useful depends entirely on the relative cost of a missed positive and a false alarm.

55. It is a variant for which the available evidence is insufficient to classify as pathogenic or benign. The classification describes the state of knowledge rather than a property of the variant, which is why reclassification is routine, why laboratories re-evaluate stored results, and why about a quarter of the TP53 records in ClinVar carry this label despite TP53 being among the most intensively studied genes in medicine.

56. Confirm the call itself — depth, allele balance, mapping quality, and whether it lies in a region short reads resolve badly. Fix the transcript and the genome build, and express the variant in standard nomenclature against them. Check population frequency in a reference matched as closely as possible to the patient’s ancestry, since a severe early-onset condition cannot be caused by a variant carried by 1% of people. Assemble consequence predictions, domain context and any functional data; consult ClinVar with attention to review status rather than label, and the literature. Then apply the ACMG/AMP criteria and consider segregation in the family. No clinical statement is made before confirmation in an accredited laboratory.

57. It quantifies PVS1, the criterion that assigns very strong evidence of pathogenicity to a null variant in a gene where loss of function is a known disease mechanism. The missense figure is not evidence that missense variants are usually benign: it is an average over a class that mixes the hotspot substitutions in the DNA-binding domain, which are almost all pathogenic, with substitutions elsewhere in the protein, which mostly are not. An average over a heterogeneous class tells you about the class, not about a variant in it.

58. It supports the claim that pathogenic missense variants in this gene concentrate in a specific structural region far beyond what its length would predict, which is the population-level evidence behind PM1 — a variant in a mutational hotspot or a critical functional domain. It does not tell you that a particular variant inside the domain is pathogenic; location is one moderate piece of evidence among several, and 26.1% of benign variants also fall inside the same domain.

59. Because the genotype is itself the identifier, so there is no field to remove. The chapter’s simulation makes the claim quantitative: every one of 10,000 simulated individuals is unique at 25 independent common markers, and the random match probability falls below one in eight billion at 30 markers. The governing concept is therefore controlled access rather than anonymisation.

60. Safe projects, safe people, safe settings, safe data and safe outputs.

61. It enables an allele-frequency membership attack. The attacker needs the published case-group frequencies, allele frequencies for a matched reference population, and the genotypes of the individual being tested; because the case group is defined by the condition, confirming membership discloses the diagnosis. Defensible modifications: publish only the variants that reach the paper’s significance threshold rather than the genome-wide table, which removes most of the accumulating evidence; or place the full summary table behind the same controlled access as the individual data and cite the accession in the paper.

62. The score’s weights were estimated for markers that tag unobserved causal variants through the linkage disequilibrium structure of the discovery population. Where that structure differs, the tags no longer track the causal variants and accuracy declines with genetic distance, so the score misclassifies risk differentially by ancestry. Deployment requires validation in a sufficiently powered cohort from the target population, reporting both discrimination and calibration, with weights re-estimated or at least re-calibrated locally.

63. Commit code, configuration, the environment specification, the workflow definition and documentation — everything a stranger needs to regenerate the results. Never commit derived results and figures, large raw inputs that are already archived under an accession, credentials, or any data under controlled access; a Git history preserves deleted files in earlier commits, so git rm does not undo the disclosure.

64. It records dependencies rather than order: which inputs and which scripts each output depends on, from which the tool computes what is stale. Listing the script itself as a prerequisite means editing the analysis invalidates its outputs, which is the dependency most often forgotten. make clean && make tests whether every result in the project can in fact be regenerated from the raw data by the code in the repository.

65. Clone the repository into an empty directory; create the environment from the pinned specification; verify the raw data against the checksums recorded in the manifest; run make clean && make; then compare the regenerated numbers and figures against those in the manuscript. A failure at any step means the project is not finished.

66. Because the comparison splits the reference genes by the one property the wrong parameter affects. Genetic code 11 reads TGA as a stop, so genes containing an internal TGA are fragmented and their true 3’ ends lost, while genes without one are unaffected — and the recovery rates, 96.6% against 52.5%, follow that split almost perfectly, which no general deficiency of the predictor would produce. The two external checks were the independent RefSeq annotation, which showed 995 predictions against 524 genes, and the coding sequence length distribution with the coding density, which showed fragmentation rather than over-prediction: a median of 336 bp against 856 bp and a density falling from 94.6% to 77.1%.

16.6 Self-assessment rubric

Scores on a question bank measure recall. What an employer, supervisor or collaborator needs to know is different: what you can be left alone with. The four levels below describe competence in terms of what you can do unsupervised, and they are cumulative — each level assumes the one before it. Assess yourself honestly against the highest level at which you could complete every item without help, and treat the first item you fail as your next piece of work.

Level 1: Supervised

You can run an analysis that someone else has specified. Concretely: you can navigate a filesystem from the command line and inspect large text files without opening them in an editor; you can recognise FASTA, FASTQ, SAM/BAM, VCF and GFF3 records and say what each field means; you can run a provided script or notebook to completion and tell whether it finished successfully; you can describe what a BLAST result or an alignment shows; and you can locate a record in GenBank, UniProt or the PDB and record its accession and version. What you cannot yet do is decide which analysis to run, or detect that a result is wrong.

Level 2: Competent

You can carry out a standard analysis on new data and notice when it has gone wrong. You can write a script of fifty to two hundred lines in Python or R that reads a biological file, computes something and writes a result; you can parameterise it rather than editing constants; you can run a short-read pipeline end to end and interpret its quality-control output; you can perform a differential expression or enrichment analysis with an established package and explain why its normalisation step exists; you can use Git to record your work and reconstruct a previous state; and you can read a tool’s documentation to find the option you need rather than copying a command from a forum. Crucially, you check results against an external expectation — a count you can predict, a positive control, a known answer — rather than trusting a clean exit code.

Level 3: Proficient

You can design an analysis for a question that nobody has specified for you, and defend the design. You can choose between methods on the basis of their assumptions rather than their popularity, and state the assumption that would invalidate each; you can diagnose a silent failure by constructing a comparison that would reveal it; you can structure a project so that make clean && make regenerates every result, with a pinned environment and a workflow definition; you can quantify uncertainty appropriately, including recognising when no legitimate confidence interval is available; you can write a methods section another group could follow and a limitations paragraph that names the specific scope of the claim; and you can judge whether a dataset’s consent and governance permit the analysis you intend. You are the person a laboratory asks when a result looks surprising.

Level 4: Independent

You can take responsibility for the computational component of a programme of work. You can decide what data to generate in the first place, trading cost against the power of the analysis it will support; you can build a pipeline that others run on data you will never see, with the validation and documentation that implies; you can review someone else’s computational work and identify the artefact they have not considered; you can supervise and teach the levels below, which requires being able to explain why a method works and not only how to run it; and you can be accountable for the ethical and legal framing of a project, including data governance, benefit sharing and the decision about what may be published. Most readers of this book will reach level 3 with practice; level 4 is a matter of years and of other people’s data depending on your judgement.

Assess by artefact, not by feeling

Levels are easier to judge from products than from confidence. Level 2 is demonstrated by a script someone else ran successfully; level 3 by a repository a stranger reproduced; level 4 by a pipeline other groups use and a review you wrote that changed someone’s analysis. If you cannot point at the artefact, you are probably one level lower than you think — which is normal, and is information rather than a verdict.

16.7 Learning pathways

Readers arrive at this book from three directions and should leave it in three different ones. The routes below name what to read again, what to build, and what to learn next; each assumes the book’s software environment and bundled data still work, which the check below confirms.

import hashlib
import pandas as pd
from pybfx import book_root

manifest = pd.read_csv(book_root() / "data" / "MANIFEST.tsv", sep="\t")
mismatched = []
for row in manifest.itertuples():
    path = book_root() / "data" / row.path
    digest = hashlib.sha256(path.read_bytes()).hexdigest() if path.exists() else "missing"
    if digest != row.sha256:
        mismatched.append(row.path)

print(f"{len(manifest)} bundled files listed; "
      f"{len(manifest) - len(mismatched)} verified against their recorded checksum")
#> 21 bundled files listed; 21 verified against their recorded checksum
print("mismatched or missing:", mismatched if mismatched else "none")
#> mismatched or missing: none

A clean result means every example in the book will reproduce on this machine. A mismatch means a file has been edited or re-downloaded, and any result derived from it is suspect — which is the same check that appears as item eight of the capstone’s self-review.

The wet-lab biologist moving to data

Your advantage is that you know what the measurements mean and where they came from; your gap is engineering. Reread the command-line chapter until the shell is uncomfortable rather than frightening, then the Python chapter, and accept that two or three weeks of deliberate practice on small scripting problems is the price of everything that follows. Build this, in order: a script that computes summary statistics over a directory of your own laboratory’s FASTA or FASTQ files; a reanalysis of a published dataset in your own field from the raw data in a public archive, where you already know what the answer should be; and a version-controlled repository for one of your own projects, with a pinned environment and a Makefile, following the capstone.

Learn next, in this order: a workflow manager, since your analyses will be rerun more often than you expect (Köster and Rahmann 2012); the statistics of your own assay in depth, which for most wet-lab readers means count-based models and multiple testing rather than general machine learning (Love, Huber, and Anders 2014); and the Bioconductor object model, because it encodes the data-integrity discipline your experiments need (Huber et al. 2015). Read Wilson et al. (2017) twice a year. The habit that will distinguish you is the one you already have from the bench — recording exactly what was done — applied now to code rather than to a notebook.

The programmer moving to biology

Your advantage is that the tooling holds no mystery; your gaps are domain semantics and an unfamiliar standard of evidence. Reread Chapters 2, 3 and 4 slowly, and resist the instinct to treat biological file formats as merely badly designed: the awkwardness usually encodes a biological fact, as one-based inclusive coordinates and six reading frames do. Your characteristic failure mode is a technically excellent analysis of the wrong quantity — normalising counts in a way that destroys the comparison, or building a classifier that predicts library protocol. Build this: a reimplementation of pairwise alignment from the recurrence, which teaches what the scoring parameters mean; an annotation comparison like the capstone’s, which teaches what an annotation is; and a reproduction of a published figure from the authors’ deposited data, which teaches the field’s standard of evidence.

Learn next: the statistics of high-dimensional biology, particularly false discovery rate control and the validation design the data science chapter insists on (Benjamini and Hochberg 1995); enough molecular biology to interrogate a result, for which an undergraduate genetics text is a better investment than another library; and the ecosystem conventions — Bioconductor, Biopython, Bioconda — so that you build on the field’s infrastructure rather than beside it (Cock et al. 2009). The most valuable thing you can offer is software engineering discipline applied to analyses that biologists trust; the fastest way to waste it is to solve a problem the field does not have.

The clinician

Your advantage is that you know what the result has to support and what an uncertain answer costs a patient; your gaps are the pipeline that produced the variant and the statistics that qualify it. Reread the clinical chapter first, then the sequencing chapter, so that you can see which pipeline decisions a clinical call rests on — depth, mapping quality, and the repeat regions where short reads cannot resolve the sequence at all. Then reread the chapter on responsible practice, because those obligations will be yours in a way they are not for a laboratory scientist. Build this: a reanalysis of the bundled TP53 ClinVar records, grouping by review status rather than by label, until the difference between an expert-panel assertion and a single-submitter one is intuitive; and an application of the ACMG/AMP criteria to a handful of variants with the evidence written out criterion by criterion (Richards et al. 2015).

Learn next: population allele frequency resources and their ancestry composition, since frequency is the evidence that most often decides a classification and its power is unequal across populations (Karczewski et al. 2020); pharmacogenomic guidelines in the areas you prescribe in, where the evidence is already actionable (Relling and Evans 2015); and enough command-line and scripting competence to interrogate a VCF yourself rather than relying on a report. You do not need to become a bioinformatician. You do need to be able to ask the laboratory a question it cannot deflect, and to recognise a result that has outrun its evidence.

16.8 Summary

• The book’s through-line is representation: every method is a decision about how to turn a biological object into data, and every interpretation is a judgement about what that representation lost.
• The chapters are sequentially dependent. Alignment needs molecular biology, pipelines need the shell, statistics needs the data structures, and the capstone needs all of them; a gap in an early chapter surfaces as confusion in a later one.
• The expensive errors are silent ones. A wrong translation table, a recycled vector, a leaky cross-validation and an unnormalised count matrix all produce plausible output and no warning, and all are caught the same way — by an independent check.
• Recall, application and interpretation are separable competences. Most readers are strongest at recall and weakest at interpretation, which is the one that practice actually rewards.
• Interpretation questions almost always resolve to the same move: ask what else could produce this result, and what evidence would distinguish the alternatives.
• Competence is best assessed by artefacts — a script someone else ran, a repository a stranger reproduced, a pipeline others depend on — rather than by self-report.
• The three reader types need different next steps: engineering practice for the wet-lab biologist, domain semantics and evidential standards for the programmer, and pipeline literacy with statistical qualification for the clinician.
• Whatever route you take, the habits of the last two chapters — pinned environments, recorded provenance, declared limitations and governance settled before the analysis — are what make the technical skills usable by anyone other than yourself.

16.9 Exercises

The first three problems are integrative: each requires methods from several chapters and cannot be solved from one. The remainder are shorter.

1. From transcript to protein, across three resources. Using the bundled TP53 GenBank record and the bundled UniProt FASTA, extract the annotated coding sequence, translate it, and test whether the product is identical to the UniProt canonical sequence. Then state what a mismatch would have implied, listing the three most likely causes in order of probability. (Chapters 2, 3 and 6.)
2. What ClinVar actually contains. Using the bundled TP53 ClinVar records, tabulate clinical significance and review status, compute the proportion of records that are uncertain or conflicting, and determine how that proportion changes when restricted to records reviewed by an expert panel. Interpret the result as advice to a clinician who has been handed a variant report. (Chapters 9, 12 and 13.)
3. From structure to clinical significance. Using the bundled p53 core domain–DNA complex, identify the protein residues within 4 Å of the DNA in each protein chain, and compare that set with the commonly cited mutational hotspot codons 175, 245, 248, 249, 273 and 282. Explain what the overlap and the non-overlap each imply about how a missense substitution can destroy a transcription factor’s function. (Chapters 10, 11 and 13.)
4. For each of Chapters 4, 9 and 12, name the single assumption whose violation would most damage the chapter’s central method, and the diagnostic that detects the violation.
5. A collaborator sends you a results table and a figure but no code. List, in order, the four things you would request, and say what each one makes checkable.
6. Place yourself on the rubric in Section 16.6 and name the first item at the next level that you cannot currently do. Write the two-week plan that would let you do it, specifying the artefact you would produce.
7. Choose any question in the bank that you answered incorrectly. Write a short explanation of why your answer was wrong that would prevent you making the same error again — specifically, name the step in your reasoning that failed rather than the fact you did not know.
8. A reviewer writes: “the authors’ pipeline is standard and the result is plausible, but I cannot determine whether it is correct.” Using the capstone’s self-check, write the response that would answer the reviewer, naming what you would add to the manuscript and the repository.

16.10 Solutions

Solutions

1. The translated coding sequence matches the UniProt canonical sequence exactly at 393 residues, which is the check that confirms the transcript, the annotation of its CDS feature and the protein record all refer to the same isoform.

from Bio import SeqIO

transcript = SeqIO.read(book_data("sequences", "TP53-NM_000546.6.gb"), "genbank")
coding_feature = next(f for f in transcript.features if f.type == "CDS")
translated = coding_feature.extract(transcript.seq).translate(to_stop=True)
canonical = SeqIO.read(book_data("proteins", "P04637.fasta"), "fasta").seq

print(f"translated {len(translated)} aa, UniProt canonical {len(canonical)} aa, "
      f"identical: {translated == canonical}")
#> translated 393 aa, UniProt canonical 393 aa, identical: True

Had they differed, the likeliest causes in order are: the transcript is a different isoform from the canonical protein, which is the usual explanation for a length difference; the coding sequence was extracted without respecting the feature’s strand or its joined location, which is the error the GenBank chapter warns about; or the records come from different releases. A difference at a single residue would instead suggest a reference allele difference, such as the common codon 72 polymorphism.

2. Uncertain and conflicting records together make up a large share of all TP53 records, and the share falls sharply among the 191 expert-panel-reviewed ones. The advice to a clinician is that the review status field, not the significance label, determines how much weight an assertion carries.

import re
import pandas as pd

records = []
with open(book_data("variants", "tp53_clinvar.vcf")) as handle:
    for line in handle:
        if line.startswith("#"):
            continue
        info = line.split("\t")[7]
        significance = re.search(r"CLNSIG=([^;]+)", info)
        review = re.search(r"CLNREVSTAT=([^;]+)", info)
        records.append({"significance": significance.group(1) if significance else "not_provided",
                        "review": review.group(1) if review else "not_provided"})

clinvar = pd.DataFrame(records)
unclear = clinvar["significance"].str.contains("Uncertain|Conflicting")
print(f"{len(clinvar):,} records; {100 * unclear.mean():.1f}% uncertain or conflicting")
#> 3,888 records; 38.3% uncertain or conflicting
expert = clinvar["review"] == "reviewed_by_expert_panel"
print(f"expert-panel reviewed: {expert.sum()} records, "
      f"{100 * unclear[expert].mean():.1f}% uncertain or conflicting")
#> expert-panel reviewed: 191 records, 18.3% uncertain or conflicting

3. Only chain B is bound in the major groove, with twelve residues within 4 Å of the DNA; chains A and C make four peripheral contacts each and are crystallographic neighbours. Of the six hotspot codons, only 248 and 273 are in the contact set. That split is the answer to the question: a hotspot substitution can abolish function either by removing a direct DNA contact, as at Arg248 and Arg273, or by destabilising the fold or the zinc site that holds the DNA-binding surface in shape, as at 175, 245, 249 and 282 — which is why the zinc-coordinating residues are inactivating despite touching no DNA at all.

import numpy as np
from Bio.PDB import PDBParser

structure = PDBParser(QUIET=True).get_structure("1tup", book_data("structures", "1tup.pdb"))
model = structure[0]
# Standard residues only: ordered waters in the DNA chains would inflate every contact set.
dna_atoms = np.array([atom.coord for chain in "EF" for residue in model[chain]
                      if residue.id[0] == " " for atom in residue])

contacts = {}
for chain_id in "ABC":
    atoms = [atom for residue in model[chain_id] if residue.id[0] == " "
             for atom in residue]
    coords = np.array([atom.coord for atom in atoms])
    separations = np.linalg.norm(coords[:, None, :] - dna_atoms[None, :, :], axis=2)
    contacts[chain_id] = sorted({atoms[i].get_parent().id[1]
                                 for i in np.where(separations.min(axis=1) < 4.0)[0]})
    print(f"chain {chain_id}: {len(contacts[chain_id])} residues within 4 A of DNA")
#> chain A: 4 residues within 4 A of DNA
#> chain B: 12 residues within 4 A of DNA
#> chain C: 4 residues within 4 A of DNA

hotspots = [175, 245, 248, 249, 273, 282]
contacting = sorted(set().union(*contacts.values()))
print("hotspot codons contacting DNA:    ", [h for h in hotspots if h in contacting])
#> hotspot codons contacting DNA:     [248, 273]
print("hotspot codons not contacting DNA:", [h for h in hotspots if h not in contacting])
#> hotspot codons not contacting DNA: [175, 245, 249, 282]

4. Chapter 4: that similarity implies homology, violated when an alignment is short or compositionally biased; the diagnostic is the E-value read together with the aligned length, since the chapter’s shuffled decoys reached 87.5% identity over short stretches. Chapter 9: that a read originates from the locus it aligns to, violated by repeats and paralogues; the diagnostics are MAPQ, coverage uniformity and, where ground truth exists, the read’s own simulated origin. Chapter 12: that the validation data had no influence on model fitting, violated by feature selection or normalisation performed before splitting; the diagnostic is label permutation, which must drive performance to chance.

5. The analysis code, which makes the computation checkable; the environment specification with versions, which makes the numbers reproducible; the raw inputs or their accessions and checksums, which makes the provenance checkable; and the parameters used at every step including defaults, which makes the decisions checkable. A fifth, worth adding if the table is a summary, is the intermediate file it was generated from.

6. No model answer. A good plan names one artefact — a repository, a script someone else runs, a reproduced published figure — with a date, rather than a quantity of reading.

7. No model answer. The instruction to name the failed reasoning step rather than the missing fact is the point: “I compared percentage identity instead of the E-value” is actionable, whereas “I did not know about E-values” is not.

8. Add to the manuscript: software versions and all parameters including defaults; an explicit statement of the metric and what it does and does not measure; a limitations paragraph naming the scope of the claim; and a data and code availability statement giving the archive accession and an archived repository version. Add to the repository: the workflow definition, the pinned environment, the input checksums, and a demonstration that a clean clone regenerates the manuscript numbers. The reviewer’s complaint is not about the method but about the absence of anything to check, and every item above converts an assertion into something checkable.

16.11 Further reading

• Huber et al. (2015) describes the Bioconductor project’s design and the object model behind SummarizedExperiment; it is the best single explanation of why the ecosystem is organised as it is.
• Wickham et al. (2019) sets out the tidyverse’s design principles, which explains why its data manipulation and plotting verbs compose the way they do.
• Cock et al. (2009) is the Biopython reference and the natural companion to further Python work on sequence and structure data.
• Pedregosa et al. (2011) documents scikit-learn, whose cross-validation and pipeline objects are the practical answer to the information-leakage problem in the data science chapter.
• Wilson et al. (2017) is the practical statement of project organisation that the capstone follows; it is short enough to reread whenever a project begins.
• Sandve et al. (2013) gives ten rules for reproducible computational research and works well as a final checklist before submission.
• Richards et al. (2015) is the ACMG/AMP variant classification standard, and is the document a clinical reader will return to most often.
• Szklarczyk et al. (2022) documents STRING, including the evidence channels whose separation the network chapter depends on.

Benjamini, Yoav, and Yosef Hochberg. 1995. “Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing.” Journal of the Royal Statistical Society Series B: Statistical Methodology 57 (1): 289–300. https://doi.org/10.1111/j.2517-6161.1995.tb02031.x.

Cock, Peter J. A., Tiago Antao, Jeffrey T. Chang, Brad A. Chapman, Cymon J. Cox, Andrew Dalke, Iddo Friedberg, et al. 2009. “Biopython: Freely Available Python Tools for Computational Molecular Biology and Bioinformatics.” Bioinformatics 25 (11): 1422–23. https://doi.org/10.1093/bioinformatics/btp163.

Huber, Wolfgang, Vincent J Carey, Robert Gentleman, Simon Anders, Marc Carlson, Benilton S Carvalho, Hector Corrada Bravo, et al. 2015. “Orchestrating High-Throughput Genomic Analysis with Bioconductor.” Nature Methods 12 (2): 115–21. https://doi.org/10.1038/nmeth.3252.

Karczewski, Konrad J., Laurent C. Francioli, Grace Tiao, Beryl B. Cummings, Jessica Alföldi, Qingbo Wang, Ryan L. Collins, et al. 2020. “The Mutational Constraint Spectrum Quantified from Variation in 141,456 Humans.” Nature 581 (7809): 434–43. https://doi.org/10.1038/s41586-020-2308-7.

Köster, Johannes, and Sven Rahmann. 2012. “Snakemake—a Scalable Bioinformatics Workflow Engine.” Bioinformatics 28 (19): 2520–22. https://doi.org/10.1093/bioinformatics/bts480.

Love, Michael I, Wolfgang Huber, and Simon Anders. 2014. “Moderated Estimation of Fold Change and Dispersion for RNA-Seq Data with DESeq2,” February. https://doi.org/10.1101/002832.

Pedregosa, Fabian, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathieu Blondel, et al. 2011. “Scikit-learn: Machine Learning in Python.” Journal of Machine Learning Research 12: 2825–30.

Relling, Mary V., and William E. Evans. 2015. “Pharmacogenomics in the Clinic.” Nature 526 (7573): 343–50. https://doi.org/10.1038/nature15817.

Richards, Sue, Nazneen Aziz, Sherri Bale, David Bick, Soma Das, Julie Gastier-Foster, Wayne W. Grody, et al. 2015. “Standards and Guidelines for the Interpretation of Sequence Variants: A Joint Consensus Recommendation of the American College of Medical Genetics and Genomics and the Association for Molecular Pathology.” Genetics in Medicine 17 (5): 405–24. https://doi.org/10.1038/gim.2015.30.

Sandve, Geir Kjetil, Anton Nekrutenko, James Taylor, and Eivind Hovig. 2013. “Ten Simple Rules for Reproducible Computational Research.” Edited by Philip E. Bourne. PLoS Computational Biology 9 (10): e1003285. https://doi.org/10.1371/journal.pcbi.1003285.

Szklarczyk, Damian, Rebecca Kirsch, Mikaela Koutrouli, Katerina Nastou, Farrokh Mehryary, Radja Hachilif, Annika L Gable, et al. 2022. “The STRING Database in 2023: Protein–Protein Association Networks and Functional Enrichment Analyses for Any Sequenced Genome of Interest.” Nucleic Acids Research 51 (D1): D638–46. https://doi.org/10.1093/nar/gkac1000.

Wickham, Hadley, Mara Averick, Jennifer Bryan, Winston Chang, Lucy McGowan, Romain François, Garrett Grolemund, et al. 2019. “Welcome to the Tidyverse.” Journal of Open Source Software 4 (43): 1686. https://doi.org/10.21105/joss.01686.

Wilson, Greg, Jennifer Bryan, Karen Cranston, Justin Kitzes, Lex Nederbragt, and Tracy K. Teal. 2017. “Good Enough Practices in Scientific Computing.” Edited by Francis Ouellette. PLOS Computational Biology 13 (6): e1005510. https://doi.org/10.1371/journal.pcbi.1005510.