Deep learning for genomics and biomedicine: sequence models, structure, and causality

Why this subject matters now

Five years ago a genomics practitioner was expected to know a convolutional network for a single regulatory task, a random forest for a variant-classification table, and a pipeline of command-line tools for alignment and calling. The center of gravity has moved. Three things changed at once. First, structure prediction was solved to useful accuracy: AlphaFold2 (Jumper et al., 2021) turned a fifty-year open problem into a forward pass, and the predicted-structure database now covers essentially every protein in UniProt. Second, unsupervised sequence models became the substrate: protein language models trained only to fill in masked residues (ESM-2, Lin et al., 2023) learn representations from which contacts, stability, and function can be read off, and the same masked-language recipe applied to DNA underlies variant scoring without a single labeled pathogenic example. Third, generative design became real: RFdiffusion (Watson et al., 2023) and ProteinMPNN (Dauparas et al., 2022) invert the structure-prediction map, designing sequences for a target backbone that fold and bind in the wet lab.

The practitioner today is expected to hold all of this in one frame. They should be able to explain why an Enformer forward pass over 200 kilobases is feasible at all when the naive attention score matrix at that length would not fit in an 80GB accelerator; why a genome-wide association study over a million variants produces thousands of genome-wide significant hits from zero causal variants if population structure is left uncorrected; why a Gaussian likelihood is the wrong noise model for single-cell counts; and why, in a field where nearly every measurement is observational, the difference between correlation and causation is not a statistical footnote but the entire scientific problem. This page builds each of those arguments from first principles.

Core theory

The central dogma, stated for engineers

The molecular biology that matters for modeling compresses to a small set of facts. DNA is a double-stranded polymer over the alphabet \(\{\mathrm{A},\mathrm{C},\mathrm{G},\mathrm{T}\}\); the two strands are complementary, A pairing with T and C with G, so one strand determines the other by the reverse complement map (read the opposite strand in the opposite direction). A gene is a stretch of DNA that is transcribed into RNA (the same alphabet with U in place of T), and a protein-coding RNA is then translated three nucleotides at a time, each codon specifying one of twenty amino acids, so a protein is a string over a twenty-letter alphabet. This is the central dogma: DNA \(\to\) RNA \(\to\) protein.

  strand (5'->3')   A  C  G  T  G  G  A  A  T  C  ...
  complement (3'->5') T  G  C  A  C  C  T  T  A  G  ...
  reverse complement  ... G  A  T  T  C  C  A  C  G  T   (read opposite strand, opposite direction)

  DNA  --transcription-->  RNA  --translation-->  protein
  ACGT... (4 letters)      ACGU...                20 amino acids
                                                  |
   regulation: which genes are transcribed, when, how much
   is set by transcription factors binding short DNA "motifs"
   in promoters and distal enhancers, sometimes 10^5 bp away.
          

The part a model spends most of its time on is regulation: only a fraction of the genome codes for protein, and the rest includes the regulatory grammar that decides which genes are expressed in which cell, at what level, at what time. Regulation is mediated by proteins called transcription factors that bind short, degenerate DNA patterns called motifs, typically six to twenty base pairs, located in promoters next to genes and in enhancers that can sit tens or hundreds of kilobases away. Two consequences drive everything downstream. Because a transcription factor reads DNA as a sequence pattern, the classical model of its binding preference is a matrix over positions, which we will show is a convolution filter. Because an enhancer can act at long range, a model of expression needs a receptive field long enough to place enhancer and promoter in the same window, which is the pressure that produced dilated convolutions and then attention.

DNA as sequence: one-hot encoding

A length-\(L\) sequence \(s \in \{\mathrm{A},\mathrm{C},\mathrm{G},\mathrm{T}\}^L\) is encoded as a binary matrix \(X \in \{0,1\}^{4\times L}\), one row per base, with \(X_{a,i}=1\) exactly when the base at position \(i\) is \(a\). This is the natural encoding because it is the one under which linear scoring is meaningful: a linear functional \(\langle W, X_{\cdot,\,i:i+w}\rangle\) reads off, position by position, the weight assigned to whichever base is actually present, and nothing else. Ordinal encoding (A=0, C=1, G=2, T=3) would smuggle in a false claim that G is "between" C and T. Reverse complement is a clean operation in one-hot space: with the alphabet ordered ACGT, complementation A\(\leftrightarrow\)T, C\(\leftrightarrow\)G is exactly reversal of the channel axis, and strand reversal is reversal of the length axis, so the reverse complement of \(X\) is \(X\) flipped along both axes. This fact reappears below as an exact equivariance the architecture can be built to respect.

Position weight matrices: the classical motif model

Before any neural network, the standard model of a transcription factor's sequence preference was the position weight matrix. Collect aligned instances of the bound site and estimate, for each column \(j \in \{0,\dots,w-1\}\), the probability \(p_{a,j}\) of base \(a\) at that position; the matrix \(P = (p_{a,j})\) is the position probability matrix. Under the standard assumption that columns are independent, the likelihood of a window \(x_{0:w}\) under the motif model relative to a background model \(b\) (base frequencies \(b_a\), often uniform at \(1/4\)) is the product of per-column odds, and taking logarithms turns the product into a sum:

The matrix \(W \in \R^{4\times w}\) of log-odds is the position weight matrix proper. A positive entry says a base is enriched over background at that column; a negative entry says it is depleted. Scanning a sequence for the motif means computing \(\text{score}\) at every start position and thresholding. The independence-of-columns assumption is the model's main limitation and the reason neural networks improved on it: real binding sites have position dependencies a rank-one-per-column additive score cannot express. But the additive score is exactly what makes the scan a convolution.

Information content: bits per position

How informative is a motif column? The natural currency is the reduction in uncertainty relative to a uniform background. A uniform distribution over four bases has entropy \(\log_2 4 = 2\) bits; an observed column with distribution \(p_{\cdot,j}\) has entropy \(H_j = -\sum_a p_{a,j}\log_2 p_{a,j}\). The information content of the column is the difference,

A fully conserved column (one base with probability one) has \(H_j=0\) and \(\mathrm{IC}_j=2\) bits; a uniform column has \(\mathrm{IC}_j=0\). This is the height of the letters in a sequence logo: the column is drawn at total height \(\mathrm{IC}_j\), each letter scaled by its probability. Total information content is a scale-free measure of how specific a motif is, and it predicts how often the motif appears by chance: a motif with \(I\) bits of information is expected roughly once every \(2^{I}\) base pairs of random sequence. The synthetic motif planted in the convolution experiment on this machine has total information content \(3.997\) bits over \(8\) columns; the stronger motif used in the motif-recovery experiment has \(9.26\) bits, so it is expected by chance about once per \(2^{9.26}\approx 613\) bp, comfortably rare inside a \(200\) bp window.

Why a PWM scan is a one-dimensional convolution

Write the scan score at start position \(s\) using the one-hot matrix. Since \(x_{s+j}\) is the base at position \(s+j\), and \(X_{a,s+j}=1\) exactly when \(a=x_{s+j}\),

The right-hand side is, term for term, the definition of a valid (no-padding) cross-correlation of the four-channel input \(X\) with a single filter of shape \(4\times w\): input channels indexed by \(a\), kernel positions indexed by \(j\), output position indexed by \(s\). Deep-learning libraries implement "convolution" as cross-correlation (they do not flip the kernel), so a PWM scan is exactly a Conv1d(in_channels=4, out_channels=1, kernel_size=w) whose single weight is the log-odds matrix \(W\) and whose bias is zero. There is nothing approximate here: the two computations produce the same number at every position, up to floating-point rounding. On this machine, scoring \(64\) random sequences of length \(512\) against an \(8\)-column PWM, the explicit triple loop and F.conv1d agree to a maximum absolute difference of \(6.65\times 10^{-6}\) and a maximum relative difference of \(1.4\times 10^{-7}\) (float32), while the convolution is \(20.4\times\) faster than the Python loop even at this trivial size. The reverse-complement scan is the same identity applied to the flipped input: convolving \(X\) with \(W\) flipped along both channel and position axes reproduces the scan of the opposite strand to \(7.6\times 10^{-6}\).

This equivalence is the conceptual hinge of the entire regulatory-genomics literature. A convolutional layer with many filters is a bank of learnable PWMs, and everything a classical motif scan does, a first convolutional layer subsumes; what the network adds is that the filters are learned discriminatively rather than counted from known sites, that a nonlinearity and pooling turn per-position scores into a presence detector, and that stacking layers builds motif combinations that a single PWM cannot. When people say a trained genomics CNN "learns motifs", this identity is why the claim is literal: the first-layer filters, read as matrices, align to known binding motifs. The motif-recovery experiment below demonstrates exactly that on planted data.

Convolutional models for regulatory genomics

The first wave of deep learning in genomics simply took the PWM-as-convolution insight seriously and stacked it. DeepBind (Alipanahi et al., 2015) trained a shallow convolutional network to predict, from sequence alone, the binding of a single transcription factor or RNA-binding protein, and showed that the learned filters recover known motifs and that the model scores the effect of a mutation by the change in its output. DeepSEA (Zhou and Troyanskaya, 2015) scaled the idea to a multi-task setting: one convolutional trunk over a 1 kb window with hundreds of output units, each predicting a chromatin feature (transcription factor binding, DNase hypersensitivity, histone marks) measured by the ENCODE and Roadmap consortia. Basset (Kelley et al., 2016) did the same for chromatin accessibility across many cell types. The architecture is always the same shape: a stem convolution acting as a learnable motif bank, then alternating convolution and pooling to build hierarchical features and shrink the length, then a dense head with one logistic output per task.

The multi-task framing is not incidental. Chromatin assays are noisy and many tasks share the same underlying motifs, so predicting hundreds of tracks jointly regularizes the trunk toward features that generalize across them, and the shared trunk is what makes a single sequence-to-function model worth training. It also reframes variant interpretation: to ask what a single-nucleotide change does, feed the reference and alternate sequences through the same trunk and read the difference in every output track, a computation the field calls in silico mutagenesis, examined in detail below.

Dilated convolutions and the receptive-field arithmetic

A plain convolutional stack grows its receptive field linearly in depth, which is hopeless for enhancer-promoter distances of \(10^4\) to \(10^5\) bp: reaching a \(10^5\) bp field with width-\(3\) convolutions would take about \(50{,}000\) layers. Two devices fix this. Pooling multiplies the spacing between output positions (the jump), so each subsequent convolution reaches farther in input coordinates. Dilated convolutions (Yu and Koltun, 2016) space their taps apart by a dilation factor \(d\), so a width-\(k\) dilated filter spans \((k-1)d+1\) input positions while using only \(k\) weights; stacking dilations that double, \(1,2,4,8,\dots\), grows the receptive field geometrically. This is the Basenji and Enformer trunk (Kelley et al., 2018; Avsec et al., 2021): a convolutional stem and a few pooling blocks to reduce the length and build local motif features, then a tower of dilated residual convolutions to reach across the whole locus.

The receptive field is computed by a two-line recurrence. Let \(r_\ell\) be the receptive field (in input base pairs) of one output unit after layer \(\ell\), and \(j_\ell\) the jump (input spacing between adjacent output units). A layer with kernel size \(k_\ell\), stride \(s_\ell\), dilation \(d_\ell\) updates them by

The jump is a product of strides because each stride-\(s\) layer keeps one of every \(s\) positions. The receptive-field increment is \((k_\ell-1)\,d_\ell\,j_{\ell-1}\) because the filter has \(k_\ell-1\) taps beyond its center, each dilated tap reaches \(d_\ell\) output-grid steps away, and one output-grid step at this depth is \(j_{\ell-1}\) input base pairs. Applying the recurrence to the trunk benchmarked on this machine (a width-\(11\) stem, four width-\(5\) convolutions each followed by \(2\times\) pooling, then a dilated tower with \(d\in\{1,2,4,8,16,32,64\}\) and kernel \(3\)) gives the trace below.

The analytic receptive field is \(8290\) bp and the output bin size (the final jump) is \(32\) bp, so each output position summarizes a \(32\) bp bin and depends on \(8290\) bp of input around it. The script confirms the arithmetic empirically by a gradient check: it feeds a real-valued input, backpropagates from one central output unit, and counts input positions with nonzero gradient. On this machine the measured support saturated at \(4096\), which is exactly the input length used in that check; the analytic field of \(8290\) bp exceeds the \(4096\) bp input, so the empirical support is clipped by the sequence boundary rather than contradicting the arithmetic. Feeding a longer input than the receptive field is required to measure the full field, a boundary subtlety worth internalizing: the receptive field a model advertises is only realized when the input is at least that long. Enformer stretches this idea to a \(196{,}608\) bp (\(\approx\!200\) kb) input with a \(128\) bp output bin, which is where the conv trunk hands off to attention.

Transformers at genomic scale, and why context length is the constraint

Dilated convolutions reach far but sum contributions with fixed, position-independent weights; they cannot, in one layer, route information adaptively from a specific distal enhancer to a specific promoter based on content. Attention can: it computes, for every pair of positions, a data-dependent weight and mixes accordingly. Enformer (Avsec et al., 2021) keeps the convolutional stem and pooling to reduce a \(196{,}608\) bp input to \(1536\) bins of \(128\) bp, then applies transformer layers over those bins, and the attention improved long-range enhancer-to-gene predictions over the purely convolutional Basenji. The reason the conv trunk stays is arithmetic: self-attention over a length-\(N\) sequence forms an \(N\times N\) score matrix, so both compute and the memory to materialize the scores scale as \(N^2\). This is the same wall analyzed in the attention and language-models pages; genomics just hits it at brutal lengths, because a genome is longer than any document.

The measured cost on this machine (an H100 80GB, bf16, \(8\) heads of dimension \(64\), batch \(1\)) makes the wall concrete. The score matrix alone, at \(H\!\cdot\!N^2\) entries in bf16, is \(0.27\) GB at \(N=4096\), \(4.3\) GB at \(N=16384\), \(68.7\) GB at \(N=65536\), and \(618\) GB at \(N=196608\) (Enformer's input length before pooling). Materializing that matrix, the "naive" path, runs out of memory beyond \(N=16384\) on an 80GB card; a fused kernel that never writes the scores (FlashAttention, Dao et al., 2022) runs the same math with peak memory that stays under \(1\) GB at every length.

Two lessons sit in this table. First, the flash kernel does not change the \(N^2\) compute (both paths do \(4HN^2D\) FLOPs), only the \(N^2\) memory; it is why a \(200\) kb attention is runnable at all, not why it is cheap. At \(N=196608\) it still takes \(255\) ms per forward attention. Second, contrast the dilated conv trunk over the same lengths: on this machine its forward pass is essentially flat, \(0.648\) ms at \(N=16384\), \(0.655\) ms at \(65536\), \(0.657\) ms at \(196608\), with peak memory under \(0.16\) GB, because convolution is linear in \(N\). That gap, quadratic attention against linear convolution, is precisely why genomic architectures put convolutions first to shrink \(N\) by pooling and reserve attention for a manageable number of bins. Context length, not parameter count, is the binding constraint in genomic modeling, and it is why so much recent work (long-range attention variants, state-space models, Hyena/StripedHyena in Evo, Nguyen et al., 2024) targets sub-quadratic sequence mixing specifically for DNA.

Variant effect prediction: scoring a mutation

The clinical payload of a sequence model is a number for a mutation: is this single-nucleotide change benign or damaging? There are two families of method. A supervised task model like DeepSEA or Enformer scores a variant by in silico mutagenesis: run the reference sequence and the alternate (mutated) sequence through the same trunk and take the difference in the predicted functional tracks, \(\Delta = f(x_{\text{alt}}) - f(x_{\text{ref}})\). A large \(\Delta\) in an accessibility or binding track says the variant disrupts a regulatory element. An unsupervised sequence model scores a variant with no labels at all, by the log-likelihood ratio the model assigns to the alternate versus reference base in its context:

The logic is evolutionary: a model trained to predict held-out residues or bases across many sequences learns which substitutions are tolerated where; a variant the model finds surprising (large negative LLR relative to reference) sits at a conserved, constrained position and is a candidate pathogenic change. This is the mechanism behind ESM-1v's zero-shot variant effect scoring for proteins (Meier et al., 2021) and AlphaMissense's pathogenicity predictions (Cheng et al., 2023), and it is evaluated by how well the LLR ranks known pathogenic variants above benign ones, typically by area under the ROC curve against ClinVar or deep mutational scanning assays.

The synthetic experiment on this machine shows both the promise and the honest limits of the LLR approach. A small three-layer masked transformer (\(d=128\), \(4\) heads, \(604{,}164\) parameters) is trained only on masked-base prediction over \(64\) bp sequences that each carry one draw from an \(8\) bp motif in random background. It reaches a masked cross-entropy of \(1.3856\) nats, barely below the uniform baseline \(\log 4 = 1.3863\) nats, because in this synthetic world almost the entire sequence is irreducibly random and only the \(8\) motif positions are predictable. Read off with no labels, the LLR is systematically more negative for substitutions inside the planted motif (mean \(-0.0876\)) than outside it (mean \(-0.0007\), standard deviation \(0.155\)), and ranking variants by \(-\mathrm{LLR}\) separates in-motif from out-of-motif changes with AUROC \(0.653\). The signal is real but weak here precisely because the context is random noise; the same recipe works far better on real genomes and proteins, where evolutionary context is rich and constrained, which is the whole point of training these models on nature's data rather than uniform sequence.

Interpretability where the science is the point: attribution and in-silico mutagenesis

In most machine learning, interpretability is a nicety; in genomics it is frequently the deliverable, because the question is which base drives the phenotype, not merely what the phenotype is. Two tools dominate. In silico mutagenesis is exhaustive: for each position, substitute each of the other three bases and record the change in the model's output, producing a \(4\times L\) matrix of effects that is read exactly like a PWM, with the sign and magnitude at each base showing what the model thinks that position contributes. Gradient-based attribution (saliency, and the axiomatic integrated gradients of Sundararajan et al., 2017, together with DeepLIFT, Shrikumar et al., 2017) approximates the same object in one backward pass instead of \(3L\) forward passes, by attributing the output to each input base.

The motif-recovery experiment ties this back to the PWM identity and closes the loop. On this machine a small convolutional classifier (\(32\) filters of width \(12\), \(1601\) parameters) is trained to detect whether a \(200\) bp sequence contains one draw from a planted \(8\) bp, \(9.26\)-bit motif, reaching test AUROC \(0.897\); a control trained on shuffled labels reaches \(0.506\), i.e. chance, confirming the signal is in the motif and not an artifact. Reading the trained filters as matrices, \(11\) of the \(32\) filters recover the planted PWM with Pearson correlation above \(0.9\) and the best filter matches at \(r=0.998\); the convolution really did learn the generating motif. And the in silico mutagenesis map averaged over positives, the model's own attribution, correlates with the true PWM log-odds at \(r=0.820\). The learned filter, the attribution map, and the classical PWM are three views of the same object, which is what the PWM-as-convolution derivation promised. For comparison, a pre-2015-style \(6\)-mer logistic-regression baseline over \(4096\) features reaches AUROC \(0.769\) on the same task: the convolution's advantage is that it detects the motif regardless of position and does not blow up its feature count with motif width.

Protein language models: masked prediction on amino-acid strings

A protein is a string over twenty letters, and the transformer that works for text works for proteins with the alphabet swapped. ESM-2 (Lin et al., 2023, at Meta AI) is a BERT-style masked language model trained on tens of millions of protein sequences from UniRef: mask a fraction of residues, predict them from the rest, scale to billions of parameters. The architecture is the encoder transformer derived in the language-models-from-scratch page, so the machinery, tokenization, positional encoding, multi-head attention, the masked cross-entropy objective, transfers directly; the change is the data distribution. What the embeddings capture is striking: without any structural supervision, the model's internal representations encode secondary structure, and, most importantly, the attention maps of specific heads light up on residue contacts, pairs of residues far apart in sequence but close in the folded structure. A linear probe on the attention maps recovers a contact map, and ESMFold (the same paper) turns the language-model embeddings into a folding head that predicts structure directly from a single sequence, trading some accuracy for the removal of the multiple-sequence-alignment search that AlphaFold requires.

Why should a model trained only to fill in masked residues learn contacts? The next subsection answers it: contacts leave a statistical fingerprint in the distribution of protein sequences, the coevolution signal, and any model that predicts a residue well from its context must internalize that fingerprint. The protein language model discovers direct-coupling structure implicitly across all proteins at once, where classical coevolution analysis estimated it per-family from an alignment.

Structure prediction I: the coevolution signal and direct coupling analysis

Two residues that touch in the folded protein are under a joint evolutionary constraint: a destabilizing mutation at one can be compensated by a mutation at the other, so across evolution the two positions covary. Given a multiple sequence alignment (MSA) of a protein family, one column per position and one row per homolog, correlated columns are evidence of contact. The naive version, raw mutual information between columns, fails because correlation is transitive: if position 1 contacts 2 and 2 contacts 3, then 1 and 3 will appear correlated through 2 even if they never touch. Direct coupling analysis (Morcos et al., 2011; Marks et al., 2011) removes these indirect chains by fitting a joint model of the whole alignment rather than scoring pairs independently. Model the probability of an alignment row \(\sigma=(\sigma_1,\dots,\sigma_N)\) as a Potts (maximum-entropy) model matching the observed single- and pairwise-column frequencies,

where \(h_i\) are per-position fields and \(J_{ij}\) are pairwise couplings. The direct coupling \(J_{ij}\) between positions \(i\) and \(j\), unlike their marginal correlation, is large only when they constrain each other after accounting for all paths through other positions; a scalar summary of \(\lVert J_{ij}\rVert\) (the direct information) ranks pairs by contact propensity. This maximum-entropy, condition-on-everything move is exactly the difference between a marginal correlation and a partial correlation, and it is the same distinction that separates a spurious association from a direct effect in the causal-inference sections below. Coevolution-derived contacts were the input that made the first deep contact predictors (RaptorX-Contact, Wang et al., 2017) and then AlphaFold work.

Structure prediction II: AlphaFold2, conceptually

AlphaFold2 (Jumper et al., 2021, at DeepMind) predicts a protein's three-dimensional structure from its sequence and an MSA of its homologs, at accuracy competitive with experiment for many targets. Four ideas carry the architecture; the paper is the reference for the full detail.

The two representations. AlphaFold2 maintains an MSA representation (rows = homologous sequences, columns = residues) and a pair representation, a matrix indexed by residue pairs \((i,j)\) that holds the model's evolving belief about their spatial relationship. The pair representation is where the coevolution signal lives inside the network: it is the learned, refined descendant of a direct-coupling contact map.

The Evoformer. The trunk is a stack of blocks (the Evoformer) that repeatedly exchange information between the two representations. Attention over the MSA lets the model read coevolutionary statistics; attention over the pair representation, structured by triangle updates that enforce a geometric consistency condition (the constraint that distances \(d_{ij}\), \(d_{jk}\), \(d_{ik}\) must be compatible, a soft triangle inequality), refines the pairwise geometry. The MSA informs the pair, and the pair informs the MSA, block after block.

The structure module. A final module turns the abstract pair and single representations into explicit \(3\)D coordinates, representing each residue as a frame (a rotation and translation) and predicting the frames directly, with an attention mechanism that operates in \(3\)D space (invariant point attention) so the output is equivariant to global rotation and translation of the protein.

Recycling. The whole network is run several times, feeding its own output structure and representations back in as input, so the model iteratively refines a prediction it could not reach in one pass. Recycling is cheap relative to depth and contributes substantially to the final accuracy. The training objective combines a structural loss on the predicted frames (the frame-aligned point error) with auxiliary losses including a masked-MSA prediction head, so the model is partly a masked language model over alignments, which is why protein language models and AlphaFold learn overlapping signals.

Structure prediction III: AlphaFold3, diffusion decoders, and design

AlphaFold3 (Abramson et al., 2024) generalized the target from single proteins to complexes, proteins with DNA, RNA, ions, and small-molecule ligands, and replaced the structure module with a diffusion decoder: the network denoises atomic coordinates from noise conditioned on the trunk representation, the same generative mechanism derived in the diffusion-and-large-vision-models page, now applied to \(3\)D coordinates rather than image pixels. A diffusion decoder is a natural fit for structure because the target is inherently multimodal (a molecule can have several valid conformations) and diffusion samples from a distribution rather than regressing a single point, and because it removes the need for the hand-built stereochemical machinery of the AlphaFold2 structure module. ESM3 (Hayes et al., 2025) took a different route, a single multimodal masked model generating over sequence, structure, and function tokens jointly.

Design is the inverse problem: not "what structure does this sequence fold to" but "what sequence folds to this structure", and, one level up, "what structure performs this function". RFdiffusion (Watson et al., 2023, at the Baker lab) is a diffusion model over protein backbones, fine-tuned from the RoseTTAFold structure predictor, that generates novel backbones satisfying constraints such as binding a given target or scaffolding a functional site; ProteinMPNN (Dauparas et al., 2022) then solves the sequence-design step, predicting an amino-acid sequence that will fold to a given backbone, as a conditional model over residues given geometry. The two compose into a design pipeline, backbone by diffusion then sequence by MPNN, whose outputs have been validated experimentally, including de novo binders and enzymes. This is the same generative-inverse structure seen throughout the diffusion page: the forward map (sequence \(\to\) structure) is learned first, and generation inverts it under constraints.

Single-cell genomics: why counts need a count likelihood

Single-cell RNA sequencing (scRNA-seq) measures, for each of thousands of cells, an integer count of transcripts per gene: a matrix of counts, cells by genes, that is large, sparse, and noisy. The counts are small integers with many zeros (a gene may be detected in a handful of cells), and total counts per cell vary by an order of magnitude for purely technical reasons (sequencing depth). Two modeling errors are tempting and wrong. Treating the counts as Gaussian ignores that they are non-negative integers whose variance grows with the mean; ignoring the depth differences confounds biology with capture efficiency.

The right noise model is the negative binomial, which is a Poisson whose rate is itself Gamma-distributed, giving the overdispersion (variance exceeding the mean) that real counts show. The synthetic count experiment on this machine makes the case quantitatively. Sampling \(4000\) cells over \(200\) genes from a negative binomial with inverse dispersion \(\theta=1.5\), the observed zero fraction is \(0.0529\), whereas a Poisson with the same means would produce only \(0.0103\) zeros; the empirical log-variance-versus-log-mean slope is \(1.93\), whereas a Poisson forces slope exactly \(1\) (variance equals mean). On held-out entries the per-entry log-likelihoods rank as expected: negative binomial \(-4.04\), a Gaussian on raw counts \(-4.45\), and Poisson a distant \(-17.63\), the last crushed because it cannot accommodate the overdispersion. Applying the field's other habit, \(\log(1+\text{CPM})\) transformation, does stabilize the variance, dropping the log-variance-log-mean slope from \(1.93\) to \(0.168\), which is why the transform is so common even though the principled route is to model the counts directly.

scVI (Lopez et al., 2018, at Berkeley) is the model that does it directly: a variational autoencoder (derived in the deep-generative-models tradition) whose decoder parameterizes a negative binomial (or zero-inflated negative binomial) per gene, with an explicit per-cell size factor that absorbs sequencing depth and batch-effect covariates that absorb technical differences between experiments. The latent variable is a low-dimensional cell embedding for clustering and visualization; because the likelihood is a proper count model and the size factor is separated out, the embedding reflects biological state rather than depth, and the batch covariate lets the model integrate datasets from different labs into one latent space. The general lesson, which recurs in multi-omics integration (combining RNA, chromatin accessibility, protein, and spatial measurements of the same cells), is that a shared latent variable with modality-appropriate likelihoods is the workhorse: totalVI and MultiVI extend scVI's structure to several modalities, and the integration is exactly a multi-view latent-variable model.

Evaluation is where genomics models go to die: homology leakage

The single most common way a genomics model reports a good number and fails in practice is a data split that leaks. Biological sequences are not independent: genes have paralogs, individuals share haplotypes, species share ancestry, so two sequences drawn into train and test can be near-identical by descent. A random split then measures memorization, not generalization. The split-leakage experiment on this machine quantifies it. Four hundred "families" of \(25\) members each are generated by mutating a founder sequence at \(5\%\) per position, with the label a motif rule corrupted by \(20\%\) noise applied at the family level (so the best achievable accuracy on novel families is capped at \(0.8\) AUROC). Within-family sequence identity is \(0.926\); between-family identity is \(0.251\). Under a random split over sequences the model reaches AUROC \(0.998\); under a group split that holds out entire families (the sequence analog of holding out whole chromosomes) it reaches \(0.524\), an inflation of \(0.474\). A nearest-neighbor lookup exposes the mechanism: it achieves accuracy \(1.0\) under the random split (every test sequence has a near-twin in train) and \(0.512\) under the family holdout. The random-split number is almost entirely memorized family identity. The remedy is to split by a biological grouping, chromosome, gene family, or sequence-identity cluster, so that train and test are genuinely disjoint in ancestry, and it is why serious genomics benchmarks report chromosome-held-out performance.

Causal inference I: GWAS and confounding by population structure

A genome-wide association study (GWAS) tests, for each of millions of common variants, whether its genotype is associated with a trait, and reports the variants passing a stringent significance threshold. Two statistical problems dominate its interpretation, and both are about false positives.

The first is multiple testing. Testing a million variants at the conventional \(p<0.05\) yields, under the complete null of no association anywhere, about \(50{,}000\) "significant" hits by chance. The null experiment on this machine confirms it exactly: a million independent null tests over \(5000\) individuals produce \(49{,}913\) hits at \(p<0.05\), against the expected \(50{,}000\). This is why the field uses the genome-wide significance threshold \(5\times 10^{-8}\), a Bonferroni-style correction for roughly a million independent common-variant tests: at that threshold the same million null tests produce zero hits (the smallest \(p\)-value observed across the million was \(1.4\times 10^{-6}\), nowhere near \(5\times 10^{-8}\)). The apparently paranoid threshold is exactly calibrated to the number of tests.

The second, subtler problem is confounding by population structure. Allele frequencies differ between ancestral populations, and if a trait also differs between those populations for environmental or cultural reasons, then every variant whose frequency differs by ancestry will associate with the trait, with no causal role whatsoever. The stratification experiment on this machine is a clean demonstration: two populations, \(6000\) individuals, \(20{,}000\) variants with ancestry-dependent frequencies (\(F_{ST}=0.08\)), zero causal variants, and a trait shifted by \(0.6\) standard deviations between populations for purely environmental reasons. Uncorrected, the scan reports \(7896\) variants passing \(5\times 10^{-8}\), all false, and the genomic-control inflation factor \(\lambda_{GC}\) (the ratio of the median test statistic to its null expectation) is \(41.6\), a catastrophic inflation where \(1.0\) means calibrated. The standard fix (Price et al., 2006) is to compute the top principal components of the genotype matrix, which capture the axes of ancestry, and include them as covariates. The leading principal component here correlates with the true population label at \(0.999\); regressing it out drops the false hits to \(0\) and the inflation factor to \(0.924\), essentially calibrated. The lesson generalizes far beyond genetics: an unmodeled common cause of both predictor and outcome manufactures associations, and the defense is to measure and adjust for it.

Causal inference II: Mendelian randomization and the instrumental-variable argument

GWAS finds associations; the harder question is whether a modifiable exposure causes a disease. Observational associations between, say, LDL cholesterol and heart disease are confounded by everything that covaries with cholesterol: diet, exercise, socioeconomics, reverse causation. A randomized trial removes confounding by assigning the exposure at random, but trials are expensive, slow, and often unethical. Mendelian randomization exploits a randomization nature already ran: at conception, alleles are dealt to offspring essentially at random with respect to the confounders of adult life. A genetic variant that raises the exposure is therefore a randomized nudge to it, an instrument.

Derive the instrumental-variable estimator. Let \(X\) be the exposure, \(Y\) the outcome, \(U\) an unmeasured confounder of both, and \(Z\) a candidate instrument (the genotype). Posit the linear structural model

where \(\beta\) is the causal effect we want. A naive regression of \(Y\) on \(X\) does not recover \(\beta\): substituting, \(Y\) depends on \(U\) both directly (through \(\lambda\)) and through \(X\), so \(\operatorname{Cov}(X,Y)/\Var(X) = \beta + \lambda\,\Cov(X,U)/\Var(X)\), biased by the confounder. An instrument earns its name by three assumptions:

(IV1) Relevance: \(\Cov(Z,X)\neq 0\), the variant actually affects the exposure (\(\gamma\neq 0\)). This is testable, and it is why MR uses variants with a genome-wide-significant effect on the exposure. (IV2) Independence (exchangeability): \(Z\) is independent of the confounder \(U\), \(\Cov(Z,U)=0\). This is what Mendelian randomization buys: alleles assigned at conception are independent of adult lifestyle and environment. (IV3) Exclusion restriction: \(Z\) affects \(Y\) only through \(X\), with no direct path, \(\Cov(Z,\varepsilon_Y)=0\) and no \(Z\to Y\) arrow. This is the assumption most easily violated, by pleiotropy (a variant influencing the outcome through a second pathway).

Under (IV2) and (IV3), compute the covariance of \(Z\) with \(Y\):

The confounder drops out because \(Z\) is independent of it, so

the Wald ratio: the effect of the variant on the outcome divided by its effect on the exposure. In two-sample MR the two slopes come from separate GWAS (\(\beta_{ZY}\) from an outcome study, \(\beta_{ZX}\) from an exposure study), which is what makes MR so widely applicable, it needs only published summary statistics. The toy simulation on this machine confirms the algebra: with true causal effect \(\beta=0.5\), a confounder inflating the association, and a valid instrument, the confounded OLS regression of \(Y\) on \(X\) returns \(1.114\) (badly biased), while the Wald ratio \(\beta_{ZY}/\beta_{ZX} = 0.2016/0.4019 = 0.5016\) recovers the true effect. The entire method rests on assumptions (IV2) and (IV3), which cannot be fully tested; MR's craft is in defending them, using multiple independent instruments, and detecting pleiotropy (MR-Egger, median-based estimators, Bowden et al., 2015) when they fail.

Why correlation versus causation is the central problem

Every method on this page ultimately serves one goal: understanding what in the genome causes a phenotype, so that an intervention (a drug, a gene edit) will work. But almost all genomic data is observational, and observational associations are confounded by ancestry, by linkage (a causal variant drags its non-causal neighbors along, so association localizes to a region, not a base), by cell-type composition, by technical batch. A predictive model can be excellent and causally useless: a sequence-to-expression model that has learned the regulatory grammar will correctly predict that mutating a motif changes expression, but a model that has merely memorized which regions are expressed will fail the moment it is asked about a sequence it has not seen. The recurring technical answer is the same idea in different clothes, condition on the common cause, whether by triangle updates and direct couplings that strip indirect correlations from a contact map, by principal components that strip ancestry from a GWAS, or by an instrument that isolates a randomized component of an exposure. A practitioner who cannot articulate the confounder in a given genomic analysis is not yet reading the results correctly.

Implementation

The code below is the runnable core of the mechanical claims, in PyTorch and JAX where an array framework is enough and in plain Python where the point is an estimator. The first block demonstrates the PWM-as-convolution identity and the reverse-complement identity, matching the maximum-difference figures reported above. Shapes are annotated; the one-hot layout is channels-first (\(4\times L\)) to match Conv1d.

import numpy as np
import torch
import torch.nn.functional as F

ALPHA = "ACGT"  # channel order; reverse-complement = flip both axes

def onehot(seq_idx):                       # (N, L) int -> (N, 4, L) float32
    N, L = seq_idx.shape
    x = np.zeros((N, 4, L), dtype=np.float32)
    n = np.arange(N)[:, None]; p = np.arange(L)[None, :]
    x[n, seq_idx, p] = 1.0
    return x

# a position weight matrix as log-odds against a uniform background
rng = np.random.default_rng(0)
w = 8
ppm = rng.dirichlet(np.full(4, 0.35), size=w).T          # (4, w) columns sum to 1
W = np.log2(ppm / 0.25).astype(np.float32)               # (4, w) log-odds

L, N = 512, 64
seq = rng.integers(0, 4, size=(N, L))                    # (N, L)
X = torch.from_numpy(onehot(seq))                        # (N, 4, L)

# explicit PWM scan: sum of W[base, j] over the window
scan = np.zeros((N, L - w + 1))
for i in range(N):
    for s in range(L - w + 1):
        scan[i, s] = sum(W[seq[i, s + j], j] for j in range(w))

# the SAME thing as a 1-D convolution: in=4, out=1, kernel=w, weight=W
wt = torch.from_numpy(W)[None]                           # (1, 4, w) = (out, in, k)
conv = F.conv1d(X, wt).numpy()[:, 0]                     # (N, L-w+1)
print("max |scan - conv1d| =", np.abs(scan - conv).max())     # ~6.7e-6 (float32)

# reverse-complement identity: RC of one-hot = flip channel AND length axes
wrc = torch.flip(wt, dims=(1, 2))                        # RC'd filter
conv_rc = F.conv1d(X, wrc).numpy()[:, 0]
Xrc = torch.from_numpy(onehot(3 - seq[:, ::-1]))         # RC of the sequences
conv_on_rc = F.conv1d(Xrc, wt).numpy()[:, 0]
print("max |RC-filter - RC-seq| =", np.abs(conv_rc[:, ::-1] - conv_on_rc).max())

import numpy as np
import jax, jax.numpy as jnp
from jax import lax

def onehot(seq_idx):                       # (N, L) int -> (N, 4, L) float32
    N, L = seq_idx.shape
    x = np.zeros((N, 4, L), dtype=np.float32)
    n = np.arange(N)[:, None]; p = np.arange(L)[None, :]
    x[n, seq_idx, p] = 1.0
    return jnp.asarray(x)

rng = np.random.default_rng(0)
w = 8
ppm = rng.dirichlet(np.full(4, 0.35), size=w).T          # (4, w)
W = jnp.asarray(np.log2(ppm / 0.25).astype(np.float32))  # (4, w) log-odds

L, N = 512, 64
seq = rng.integers(0, 4, size=(N, L))
X = onehot(seq)                                          # (N, 4, L)

scan = np.zeros((N, L - w + 1))
Wn = np.asarray(W)
for i in range(N):
    for s in range(L - w + 1):
        scan[i, s] = sum(Wn[seq[i, s + j], j] for j in range(w))

# conv_general_dilated is cross-correlation, exactly a PWM scan
wt = W[None]                                             # (1, 4, w) = (O, I, k)
conv = lax.conv_general_dilated(
    X, wt, window_strides=(1,), padding="VALID",
    dimension_numbers=("NCH", "OIH", "NCH"))[:, 0]       # (N, L-w+1)
print("max |scan - conv| =", float(jnp.abs(jnp.asarray(scan) - conv).max()))

# reverse complement = flip both filter axes
wrc = jnp.flip(wt, axis=(1, 2))
conv_rc = lax.conv_general_dilated(X, wrc, (1,), "VALID",
    dimension_numbers=("NCH", "OIH", "NCH"))[:, 0]
Xrc = onehot(3 - seq[:, ::-1])
conv_on_rc = lax.conv_general_dilated(Xrc, wt, (1,), "VALID",
    dimension_numbers=("NCH", "OIH", "NCH"))[:, 0]
print("max |RC-filter - RC-seq| =",
      float(jnp.abs(conv_rc[:, ::-1] - conv_on_rc).max()))

The next block builds the dilated trunk whose receptive field the table derived, and measures the field empirically by backpropagating from one output unit and counting input positions with nonzero gradient. Run on an input at least as long as the field to see the full \(8290\); on a shorter input the measured support is clipped at the input length, as noted above. The Torch and JAX trunks share weights so their outputs match to \(1.5\times 10^{-5}\), the parity check the script reports.

import torch, torch.nn as nn, torch.nn.functional as F

def receptive_field(layers):
    # layers: list of (kernel, stride, dilation); returns (rf_bp, output_jump)
    r, j = 1, 1
    for k, s, d in layers:
        r = r + (k - 1) * d * j       # taps beyond center reach d*j input bp each
        j = j * s                     # jump is the running product of strides
    return r, j

# stem k=11, then 4x (conv k=5, pool /2), then dilated tower d=1..64
layers = [(11, 1, 1), (2, 2, 1)]
for _ in range(4):
    layers += [(5, 1, 1), (2, 2, 1)]
layers += [(3, 1, d) for d in (1, 2, 4, 8, 16, 32, 64)]
print(receptive_field(layers))        # (8290, 32)

class DilatedTrunk(nn.Module):
    def __init__(self, ch=16, dilations=(1, 2, 4, 8, 16, 32, 64)):
        super().__init__()
        self.stem   = nn.Conv1d(4, ch, 11, padding=5)
        self.blocks = nn.ModuleList([nn.Conv1d(ch, ch, 5, padding=2) for _ in range(4)])
        self.dil    = nn.ModuleList([nn.Conv1d(ch, ch, 3, padding=d, dilation=d)
                                     for d in dilations])
        self.head   = nn.Conv1d(ch, 1, 1)
    def forward(self, x):             # x: (B, 4, L)
        h = F.max_pool1d(F.relu(self.stem(x)), 2)
        for b in self.blocks:
            h = F.max_pool1d(F.relu(b(h)), 2)
        for d in self.dil:
            h = h + F.relu(d(h))      # dilated residual tower
        return self.head(h)           # (B, 1, L/32)

# empirical receptive field via gradient support (use L >= 8290 for the full field)
m = DilatedTrunk().double()
L = 16384
x = torch.rand(1, 4, L, dtype=torch.float64, requires_grad=True)
out = m(x); c = out.shape[-1] // 2
out[0, 0, c].backward()
g = x.grad.abs().sum(1)[0]            # (L,) sensitivity of one output unit
nz = torch.nonzero(g > 1e-14)[:, 0]
print("measured RF =", int(nz[-1] - nz[0] + 1))   # 8290 when L exceeds the field

import numpy as np
import jax, jax.numpy as jnp
from jax import lax, grad

def receptive_field(layers):
    r, j = 1, 1
    for k, s, d in layers:
        r = r + (k - 1) * d * j
        j = j * s
    return r, j

layers = [(11, 1, 1), (2, 2, 1)]
for _ in range(4):
    layers += [(5, 1, 1), (2, 2, 1)]
layers += [(3, 1, d) for d in (1, 2, 4, 8, 16, 32, 64)]
print(receptive_field(layers))        # (8290, 32)

def conv(h, w, b, pad, dil):          # NCH input, OIH weights
    y = lax.conv_general_dilated(h, w, (1,), [(pad, pad)], rhs_dilation=(dil,),
                                 dimension_numbers=("NCH", "OIH", "NCH"))
    return y + b[None, :, None]

def trunk(params, x):                 # params from a matching Torch state_dict
    h = jax.nn.relu(conv(x, params["stem.w"], params["stem.b"], 5, 1))
    h = lax.reduce_window(h, -jnp.inf, lax.max, (1, 1, 2), (1, 1, 2), "VALID")
    for i in range(4):
        h = jax.nn.relu(conv(h, params[f"b{i}.w"], params[f"b{i}.b"], 2, 1))
        h = lax.reduce_window(h, -jnp.inf, lax.max, (1, 1, 2), (1, 1, 2), "VALID")
    for i, d in enumerate((1, 2, 4, 8, 16, 32, 64)):
        h = h + jax.nn.relu(conv(h, params[f"d{i}.w"], params[f"d{i}.b"], d, d))
    return conv(h, params["head.w"], params["head.b"], 0, 1)

# empirical receptive field: gradient of one central output unit w.r.t. the input
def central(params, x):
    out = trunk(params, x); c = out.shape[-1] // 2
    return out[0, 0, c]

# g = grad(central, argnums=1)(params, x); count |g| > 1e-14 positions -> 8290

The last two blocks are pure Python (NumPy), because the point is an estimator, not a network. The first computes information content, the currency of motif specificity; the second is the Mendelian-randomization Wald estimator on a simulated confounded system, reproducing the recovery of the true causal effect \(\beta=0.5\) that the confounded regression misses.

import numpy as np

def info_content(ppm):                     # ppm: (4, w), columns sum to 1
    p = np.clip(ppm, 1e-12, 1.0)
    ic_per_col = 2.0 + np.sum(p * np.log2(p), axis=0)   # 2 - H_j, bits
    return ic_per_col, float(ic_per_col.sum())

col0 = np.array([0.7, 0.1, 0.1, 0.1])      # from Problem 1
col1 = np.array([0.85, 0.05, 0.05, 0.05])
ic, total = info_content(np.stack([col0, col1], axis=1))
print(np.round(ic, 3), "total", round(total, 3))   # [0.643 1.152] total 1.795
# expected chance spacing of a motif with I bits: about one per 2**I bp
print("chance spacing (bp):", round(2 ** total, 1))

import numpy as np
rng = np.random.default_rng(0)
n = 200_000

Z = rng.binomial(2, 0.3, n).astype(float)          # instrument: SNP dosage 0/1/2
U = rng.normal(size=n)                             # unmeasured confounder
beta = 0.5                                          # TRUE causal effect of X on Y
X = 0.4 * Z + 1.2 * U + rng.normal(0, 0.5, n)      # exposure: Z relevant, U confounds
Y = beta * X + 0.9 * U + rng.normal(0, 0.5, n)     # outcome: U confounds; Z only via X

# confounded observational estimate: regress Y on X
b_obs = np.cov(X, Y)[0, 1] / np.var(X)
# instrumental-variable (Wald ratio) estimate
b_zx = np.cov(Z, X)[0, 1] / np.var(Z)              # first stage
b_zy = np.cov(Z, Y)[0, 1] / np.var(Z)              # reduced form
b_iv = b_zy / b_zx                                  # = cov(Z,Y)/cov(Z,X)

print("true beta      :", beta)
print("OLS (confounded):", round(b_obs, 4))        # ~1.114, badly biased upward
print("Wald IV        :", round(b_iv, 4))          # ~0.5016, recovers the truth

How it is done in practice

The gap between these derivations and a deployed system is mostly data engineering and scale. A production regulatory-genomics model like Enformer trains on thousands of genomic tracks (CAGE, DNase, ChIP-seq, from ENCODE, Roadmap, and FANTOM) over the whole genome, with the split by chromosome so that the homology leakage analyzed above cannot occur, and with reverse-complement and small-shift augmentation baked in so the model respects the strand symmetry the one-hot encoding makes exact. The input is \(196{,}608\) bp, and the architecture is convolution-then-attention precisely because, as the attention cost table showed, attention over that raw length would need hundreds of gigabytes for the score matrix on this machine's H100; pooling to \(1536\) bins first is what makes it fit.

Structure prediction in practice is dominated by the multiple-sequence-alignment search, not the network: AlphaFold2's accuracy depends on finding deep, diverse homolog alignments, and the alignment search over large sequence databases can take longer than the forward pass. ESMFold's contribution was removing that search by folding from a single sequence through a protein language model, trading accuracy on shallow-MSA targets for a large speedup that made metagenomic-scale folding feasible. The AlphaFold Protein Structure Database and ESM Metagenomic Atlas together released hundreds of millions of predicted structures, changing what "look up the structure" means for a working biologist.

Single-cell pipelines in practice are a preprocessing gauntlet, quality filtering, doublet removal, normalization, before any model, and the batch integration that scVI performs is often the whole game when combining atlases from different labs. The count-likelihood point from the theory section is not academic: using a Gaussian on \(\log(1+\text{CPM})\) data is the field's common shortcut and it works for visualization, but for differential expression and for integration the negative-binomial likelihood with an explicit size factor is what keeps depth from masquerading as biology. GWAS in practice runs mixed-model methods (BOLT-LMM, SAIGE, regenie) that correct for both population structure and relatedness via a genetic random effect rather than a handful of principal components, but the principle is identical to the PCA correction demonstrated here: model the structure and condition it out, or report inflated nonsense.

The current research frontier

The frontier is being pushed on several fronts by different groups. On long-range DNA modeling, the pressure is to escape the quadratic attention wall: HyenaDNA and the Evo models (Nguyen et al., 2024, at the Arc Institute) use implicit long convolutions and state-space mixing to reach hundreds of kilobases to the megabase at single-nucleotide resolution, and the Nucleotide Transformer (Dalla-Torre et al., 2024, InstaDeep) scaled masked DNA language models across many genomes; this is the DNA analog of the state-space work covered in the sequence-models pages. On structure, AlphaFold3 (DeepMind) and its open reimplementations, along with Boltz-1 (MIT, 2024) and Chai-1, extended diffusion-based prediction to complexes and ligands, while ESM3 (EvolutionaryScale, 2025) unified sequence, structure, and function in a single generative model. On design, RFdiffusion and its all-atom successors (Baker lab, Washington), together with ProteinMPNN, moved from proof-of-concept to experimentally validated binders and enzymes, and generative design of antibodies and enzymes is the active edge. On single-cell foundation models, Geneformer (Theodoris et al., 2023, Broad) and scGPT (Cui et al., 2024, Toronto) pretrain transformers on tens of millions of cells to transfer across tasks, though whether these outperform careful task-specific models remains contested, an honest open question in the field. Across all of these the shared trajectory is the one this page traced: from supervised task models, to unsupervised sequence models whose representations transfer, to generative models that invert the forward map.

Open source to read

• facebookresearch/esm — ESM-2 and ESMFold. Open esm/pretrained.py to see the masked-language-model checkpoints, then esm/model/esm2.py for the encoder; the contact-prediction head shows how coevolution is read from attention.
• google-deepmind/alphafold and google-deepmind/alphafold3 — the reference implementations. In AlphaFold2 read alphafold/model/modules.py for the Evoformer and triangle updates; AlphaFold3 is where the diffusion structure decoder lives.
• RosettaCommons/RFdiffusion — backbone generation by diffusion. The inference config exposes the constraint types (binder, motif scaffolding) that make it a design tool rather than a sampler.
• dauparas/ProteinMPNN — sequence design for a fixed backbone. protein_mpnn_run.py is the entry point; the model is a compact message-passing network over residue geometry.
• scverse/scvi-tools — scVI and its multi-modal relatives. Read scvi/module/_vae.py for the negative-binomial decoder and the size-factor and batch handling that the count-model section derived.
• calico/basenji — the dilated-convolution regulatory-genomics trunk and the Enformer lineage. The model files make the pooling-then-dilation receptive-field construction concrete.
• instadeepai/nucleotide-transformer — masked DNA language models across many genomes, with the tokenization and probing code for downstream regulatory tasks.
• kundajelab/deeplift and shap/shap — attribution methods for reading which bases drive a prediction, the tooling behind the interpretability section.

Common misconceptions

"A convolutional genomics model learns something fundamentally different from a PWM." Its first layer is a bank of PWMs: the derivation showed the scan and the convolution are the same computation, and the motif-recovery experiment showed the trained filters align to the planted PWM at \(r=0.998\). What the network adds is discriminative fitting, position invariance through pooling, and motif combinations across layers, not a different first-layer primitive.

"FlashAttention makes long-context attention cheap." It makes it fit, not cheap. Both the naive and fused paths do the same \(O(N^2)\) floating-point work; the fused kernel only avoids materializing the \(O(N^2)\) score matrix in memory. On this machine attention at \(N=196608\) still takes \(255\) ms per forward pass even fused, which is why genomic models pool the length down before attending.

"A high test accuracy means the model generalizes." Not if the split leaks homology. The same task on this machine scored AUROC \(0.998\) under a random split and \(0.524\) under a family holdout; a nearest-neighbor lookup got \(1.0\) under the random split by pure memorization. Genomics numbers are meaningful only with a biology-aware split.

"A GWAS hit is a causal variant." It is an associated region. Linkage disequilibrium drags non-causal neighbors along with a causal variant, so association localizes to a haplotype block, not a base, and uncorrected population structure can manufacture thousands of hits with no causal variant anywhere (\(7896\) false hits in the stratification experiment). Fine-mapping and orthogonal evidence are needed to name the causal base.

"Mendelian randomization proves causation." It licenses a causal estimate only under assumptions it cannot fully test, in particular the exclusion restriction (no pleiotropy). A variant that affects the outcome through a second pathway violates it and biases the Wald ratio, as Problem 6 showed. MR is a strong argument, not a proof, and its credibility rests on defending its instruments.

"AlphaFold predicts the physics of folding." It predicts the structure, largely by exploiting the coevolutionary statistics in a multiple sequence alignment, not by simulating a folding trajectory. Its accuracy degrades when few homologs exist, which is exactly the signature of a method that leans on evolutionary covariation rather than on a force field.

"Single-cell counts can be treated as Gaussian after log transformation." The transform stabilizes variance for visualization, but the data are overdispersed counts: on this machine the observed zero fraction was five times what a Poisson predicts and the variance grew with the mean at slope \(1.93\), not \(1\). Differential expression and dataset integration want a negative-binomial likelihood with an explicit depth factor.

"Bigger sequence models are strictly better for genomics." Context length, not parameter count, is usually the binding constraint, and the benefit of a foundation model over a careful task-specific model is contested in both regulatory genomics and single-cell. The field's honest position is that scale helps unevenly and that the right architecture for the sequence length often matters more than raw size.

Self-check

References

1. Alberts, Johnson, Lewis, Morgan, Raff, Roberts, Walter. Molecular Biology of the Cell, 6th ed. Garland Science, 2014. The standard grounding for the central dogma, transcription, and regulation.
2. Lodish, Berk, Kaiser, Krieger, et al. Molecular Cell Biology, 8th ed. W. H. Freeman, 2016. Complementary textbook treatment of gene expression and its control.
3. Zou, Huss, Abid, Mohammadi, Torkamani, Telenti. A primer on deep learning in genomics. Nature Genetics, 2019. doi:10.1038/s41588-018-0295-5
4. Alipanahi, Delong, Weirauch, Frey. Predicting the sequence specificities of DNA- and RNA-binding proteins by deep learning (DeepBind). Nature Biotechnology, 2015. doi:10.1038/nbt.3300
5. Zhou, Troyanskaya. Predicting effects of noncoding variants with deep learning-based sequence model (DeepSEA). Nature Methods, 2015. doi:10.1038/nmeth.3547
6. Kelley, Snoek, Rinn. Basset: learning the regulatory code of the accessible genome with deep convolutional neural networks. Genome Research, 2016. doi:10.1101/gr.200535.115
7. Kelley, Reshef, Bileschi, Belanger, McLean, Snoek. Sequential regulatory activity prediction across chromosomes with convolutional neural networks (Basenji). Genome Research, 2018. doi:10.1101/gr.227819.117
8. Avsec, Agarwal, Visentin, Ledsam, Grabska-Barwinska, Taylor, Assael, Jumper, Kohli, Kelley. Effective gene expression prediction from sequence by integrating long-range interactions (Enformer). Nature Methods, 2021. doi:10.1038/s41592-021-01252-x
9. Yu, Koltun. Multi-scale context aggregation by dilated convolutions. ICLR, 2016. arXiv:1511.07122
10. Dao, Fu, Ermon, Rudra, Re. FlashAttention: fast and memory-efficient exact attention with IO-awareness. NeurIPS, 2022. arXiv:2205.14135
11. Nguyen, Poli, et al. Sequence modeling and design from molecular to genome scale with Evo. Science, 2024. doi:10.1126/science.ado9336
12. Dalla-Torre et al. The Nucleotide Transformer: building and evaluating robust foundation models for human genomics. Nature Methods, 2024. doi:10.1038/s41592-024-02523-z
13. Meier, Rao, Verkuil, Liu, Sercu, Rives. Language models enable zero-shot prediction of the effects of mutations on protein function (ESM-1v). NeurIPS, 2021. bioRxiv:2021.07.09.450648
14. Cheng, Novati, Pan, et al. Accurate proteome-wide missense variant effect prediction with AlphaMissense. Science, 2023. doi:10.1126/science.adg7492
15. Lin, Akin, Rao, Hie, et al. Evolutionary-scale prediction of atomic-level protein structure with a language model (ESM-2/ESMFold). Science, 2023. doi:10.1126/science.ade2574
16. Morcos, Pagnani, Lunt, et al. Direct-coupling analysis of residue coevolution captures native contacts across many protein families. PNAS, 2011. doi:10.1073/pnas.1111471108
17. Jumper, Evans, Pritzel, et al. Highly accurate protein structure prediction with AlphaFold. Nature, 2021. doi:10.1038/s41586-021-03819-2
18. Abramson, Adler, Dunger, et al. Accurate structure prediction of biomolecular interactions with AlphaFold3. Nature, 2024. doi:10.1038/s41586-024-07487-w
19. Watson, Juergens, Bennett, et al. De novo design of protein structure and function with RFdiffusion. Nature, 2023. doi:10.1038/s41586-023-06415-8
20. Dauparas, Anishchenko, Bennett, et al. Robust deep learning-based protein sequence design using ProteinMPNN. Science, 2022. doi:10.1126/science.add2187
21. Lopez, Regier, Cole, Jordan, Yosef. Deep generative modeling for single-cell transcriptomics (scVI). Nature Methods, 2018. doi:10.1038/s41592-018-0229-2
22. Price, Patterson, Plenge, Weinblatt, Shadick, Reich. Principal components analysis corrects for stratification in genome-wide association studies (EIGENSTRAT). Nature Genetics, 2006. doi:10.1038/ng1847
23. Davey Smith, Hemani. Mendelian randomization: genetic anchors for causal inference in epidemiological studies. Human Molecular Genetics, 2014. doi:10.1093/hmg/ddu328
24. Bowden, Davey Smith, Burgess. Mendelian randomization with invalid instruments: effect estimation and bias detection through Egger regression. International Journal of Epidemiology, 2015. doi:10.1093/ije/dyv080
25. Sundararajan, Taly, Yan. Axiomatic attribution for deep networks (integrated gradients). ICML, 2017. arXiv:1703.01365
26. Shrikumar, Greenside, Kundaje. Learning important features through propagating activation differences (DeepLIFT). ICML, 2017. arXiv:1704.02685