Chapter 24 — Reading BulkFormer at Author Depth

Third of the four model readings, and the one closest to home. scGPT and Tahoe both model single cells; classic toxicogenomics — DrugMatrix, Open TG-GATEs, LINCS L1000 — is almost entirely bulk: one expression profile per tissue sample, averaged over millions of cells. BulkFormer is the first foundation model built for that dense bulk modality, and it independently adopts several of the fixes the scGPT and Tahoe chapters proposed — it keeps expression magnitude (continuous values, not binning), it feeds the whole ~20,000-gene profile (no highly-variable-gene subsetting), and it initializes gene identities from protein sequence (the cross-species trick from 22_READING_SCGPT.md §22.8-4). So this chapter is partly “how a bulk FM works” and partly “how close is the field’s best bulk model to something toxicogenomics could actually use?” The honest answer, developed across §24.7–24.8: structurally very close, biologically still off on dose, time, species, and fold-change.

One-sentence version: BulkFormer is a ~147M-parameter encoder that takes an entire dense bulk transcriptome (all 20,010 human protein-coding genes, as continuous log-TPM values), mixes a gene co-expression graph with linear-attention transformer layers, and is pretrained to impute masked genes — a magnitude-preserving, whole-genome bulk model that is the right modality for toxicogenomics but trained on the wrong (observational, human, absolute-expression) data.

24.1 What problem BulkFormer solves

Bulk vs single-cell, for a non-biologist. A single-cell RNA-seq profile (Chapters 22–23) measures one individual cell: sparse, mostly zeros, ~3,000 genes detected. A bulk RNA-seq profile measures a whole tissue sample — a biopsy, a treated cell culture, an organ slice — by grinding up millions of cells and averaging their mRNA together. The result is one dense vector in which essentially every gene is nonzero (BulkFormer requires ≥14,000 detected genes per sample) and the value is a population average. Bulk is cheaper, older, and by far the dominant modality in clinical, pharmaceutical, and toxicological transcriptomics — you dose a rat, harvest its liver, and sequence the bulk.

The gap. Every prior transcriptomic foundation model — Geneformer, scGPT, scFoundation — was pretrained only on sparse single-cell data. BulkFormer’s motivating claim is that this is a modality mismatch: a model tuned to ~3,000-gene sparse cells underperforms on ~16,000-gene dense bulk samples. BulkFormer is the first FM purpose-built for bulk.

What one training example is. One bulk sample = a dense continuous vector over 20,010 protein-coding genes, values in log(TPM + 1) (TPM = “transcripts per million,” a depth-normalized expression unit). Notably — and this becomes a limitation — a training example is expression only: the authors deliberately ignore sample metadata (tissue, disease, sex, age). No condition tokens.

Input → output. In: a dense bulk profile (optionally with some genes masked). Out: (i) a sample embedding and per-gene contextual embeddings for downstream use, and (ii) imputed/denoised expression values for masked or missing genes. The pretraining task is imputation; everything else is a downstream head.

Downstream tasks it supports (all in §24.5): transcriptome imputation, disease annotation, cancer-subtype classification, patient prognosis/survival, compound-perturbation response (LINCS), drug-response/IC50 prediction (GDSC), and gene-essentiality prediction (DepMap).

24.2 How the input is represented

This is where BulkFormer quietly fixes two of scGPT’s problems, so it is worth going slowly. Start from a small bulk sample:

Two things are already different from a single cell: the vector is dense (no dropout zeros), and the values are real magnitudes, not a sparse bag.

The three embeddings per gene

Each of the 20,010 gene positions is turned into a vector by combining three pieces (compare scGPT’s three summed embeddings, §22.2, but every piece here is upgraded):

1. Gene-identity embedding — from protein sequence (ESM2). Instead of the random-then-learned gene-ID row that scGPT uses (22_READING_SCGPT.md §22.2), each gene’s identity vector is warm-started — initialized from ESM2, a protein language model, run on that gene’s canonical protein sequence, then trained further. So the gene’s representation starts out encoding what protein it makes rather than starting as noise. This is exactly the cross-species idea the scGPT chapter proposed (§22.8-4): because orthologous proteins have nearly identical sequences, this identity embedding is, in principle, species-transferable — even though BulkFormer itself only trains on human.
2. Expression-value embedding — continuous, via Rotary Expression Embedding (REE). The log(TPM+1) value enters as a continuous number, encoded with a rotary scheme (REE, an expression-space analogue of rotary positional encoding). Crucially, there is no binning — the magnitude is preserved. This directly avoids scGPT’s “51 rank bins throw away how much” problem (§22.2, §22.7).
3. Sample-context embedding. An MLP compresses the entire expression vector into one sample-level vector, broadcast across all positions, so every gene “knows” the global state of the sample it sits in.

Sequence, masking, padding, metadata

• Ordering / sequence. All 20,010 genes are one sequence — the whole transcriptome, no highly-variable-gene subsetting, no chunking. (This is only affordable because of linear attention, §24.3.) As in scGPT, gene identity replaces position; there is no positional embedding over genes.
• Masking. Pretraining masks ~15% of the gene values, replacing each with a placeholder value of −10 (a value far outside the real log-TPM range, so the model can recognize “this one is hidden”). The model predicts the masked genes’ true values. This is masked-value regression, the bulk analogue of scGPT’s GEP.
• Padding. Because every sample uses the same fixed 20,010-gene vocabulary in the same order, there is effectively no variable-length padding — a fixed dense vector every time. (Contrast scGPT’s variable ~1–3k gene subsets.)
• Metadata. None enters the model. Tissue, disease, sex, age, and — critically for us — drug, dose, and time are not represented. A training example is a bare expression vector. The authors list this as an explicit limitation, and it is the single biggest gap for toxicogenomics (§24.7–24.8).

What biology is preserved or lost

The left column is why BulkFormer is the most toxicogenomics-shaped of the four models: bulk modality, magnitude kept, whole-genome. The right column is why it is still not a toxicogenomics model: no fold-change, no dose/time, human only.

24.3 How the architecture works

BulkFormer’s distinctive move is a hybrid graph + linear-attention transformer, chosen so it can process all ~20,000 genes at once. We walk one sample through.

Why not a plain transformer

Vanilla self-attention costs $O(n^2)$ in sequence length. At $n\approx 20,010$ genes, an attention matrix is $20,010\times 20,010\approx 4\times 10^8$ entries per sample per layer — infeasible. Single-cell FMs dodge this by subsetting to ~1–3k highly-variable genes; BulkFormer instead keeps all genes and makes attention itself cheaper.

The two ingredients

Performer (linear attention). Performer approximates softmax attention with a kernel trick that makes cost scale linearly in sequence length, $O(n)$ instead of $O(n^2)$. This is the key enabler of the whole-transcriptome context. (Conceptually: instead of forming the full $n\times n$ score matrix, it factorizes attention through random feature maps.)

A gene co-expression graph + GCN. In parallel, BulkFormer builds a fixed gene–gene graph: compute the absolute Pearson correlation between every gene pair across the training corpus, discard edges below |r| = 0.4, and keep each gene’s top-20 strongest neighbors. A graph convolution (GCN) layer then lets each gene mix information with its co-expressed neighbors. This injects an explicit prior — “genes that co-vary are related” — that pure attention would have to learn from scratch.

The block. BulkFormer stacks “blocks,” each = 1 GCN layer + several Performer layers. So every block does both: local structure from the co-expression graph, global mixing from linear attention.

Sizes and shapes

Five released sizes; the main “~150M” model is the 147M variant:

End-to-end shapes for one sample ($G=20,010$, $d=640$):

A worked micro-example

Suppose CYP1A2 (a liver drug-metabolism gene) is masked (set to −10) in a liver sample. Two channels help predict it: the GCN lets it borrow from its top-20 co-expressed neighbors — other liver-enriched metabolic genes that are observed — while Performer attention lets it consult the global sample state (the context embedding says “this is a liver-like profile”). The linear head then outputs a continuous predicted log-TPM, e.g. 7.3 against a true 7.4. Biologically: the model reconstructs a gene from the company it keeps. Statistically: this is sophisticated co-expression imputation — which, as §24.7 presses, is also the ceiling of what the objective can teach.

24.4 How the model is trained

The objective

A single objective: masked-value regression. Mask ~15% of a sample’s gene values (→ −10), predict them, minimize MSE on the masked positions:

$$\mathcal{L}=\frac{1}{|\mathcal{M}|}\sum _{i\in \mathcal{M}}({\hat{x}}_i-x_i{)}^2.$$

No contrastive term, no auxiliary loss — just imputation. The intent: to reconstruct a masked gene from the rest of a dense profile, the model must internalize genome-wide co-expression structure.

The data, and what it should teach

Pretraining uses 522,769 human bulk samples assembled from GEO and ARCHS4 (public expression repositories), spanning nine physiological systems, healthy and diseased. Two corrections worth stating plainly, because they matter for this book:

• It is not GTEx/TCGA. Those tissue atlases appear only in downstream evaluation; pretraining is GEO+ARCHS4.
• It is observational. These are tissue/disease snapshots, not controlled perturbations. So, exactly as with scGPT (§22.4), the objective can teach association (co-expression) but not intervention — a point that returns for toxicogenomics.

Shortcuts and weaknesses in the objective

Shortcut 1 — imputation is “smart averaging.” Because bulk genes are strongly co-expressed, a masked gene is often well-predicted by a linear combination of its neighbors. A model can score very high on imputation MSE (BulkFormer reports PCC 0.954) while mostly learning conditional co-expression, not biology. The high score is real but easy to over-read: the baseline (co-expression regression) is already strong. Fix / honest framing: report imputation against a strong co-expression/linear baseline, and weight the evaluation toward genes that are hard to predict from neighbors.

Shortcut 2 — observational data, association only. Nothing in the objective is interventional, so the model cannot learn dose- or drug-driven causal response from pretraining; any such ability must be injected downstream (LINCS/GDSC fine-tuning). Fix: interventional pretraining (§24.8).

Shortcut 3 — the co-expression graph can bake in confounds. The GCN prior is built from corpus-wide correlations, which include batch and technical covariation, not just biology. The model can lean on edges that reflect study/platform structure. Fix: build the graph from batch-corrected data, or learn/prune edges rather than fixing them by raw Pearson.

24.5 How it is used for downstream tasks

Same pattern as before — task → concrete input → adapt → concrete output → example — condensed across the seven reported tasks. Adaptation is a task head or fine-tune on the embeddings; imputation is native.

Transcriptome imputation (the native task)

• In: a bulk profile with missing/low-quality genes.
• Out: filled-in continuous values (per-gene $\hat{x}$).
• Example / result: test PCC 0.954; external TCGA PCC 0.914. Use: rescue genes lost to low sequencing depth; harmonize datasets.

Disease annotation & cancer subtype (classification)

• In: a patient’s bulk tumor/tissue profile → sample embedding.
• Adapt: a classifier head on the embedding.
• Out: a disease/subtype label. Example: 23-disease annotation weighted-F1 0.939 (vs scGPT 0.885); 33 cancer subtypes 0.833 (vs scGPT 0.830 — a thin margin, worth noting).

Patient prognosis / survival

• In: tumor bulk profile.
• Out: risk score → survival prediction. Example: across ~10,000 patients / 33 cancers, AUROC 0.747 (modest; the authors concede “considerable room for improvement” given bulk noise). A surfaced biomarker: high H4C1 in kidney cancer ↔ 5.2× mortality.

Compound perturbation (LINCS) — the toxicology-adjacent task

• In: a baseline profile + a compound (LINCS L1000-style).
• Out: predicted perturbed profile. Example: PCC 0.493 / SCC 0.430 (beats PRnet 0.408 — but note the absolute value is modest: predicting perturbation responses remains hard, echoing 18 §18.4).

Drug-response (GDSC) & gene essentiality (DepMap)

• Drug response: cell-line bulk profile + drug → IC50; PCC 0.910 (vs scFoundation 0.880). Essentiality: predict DepMap gene-dependency; PCC 0.931.

Toxicity relevance (why this book cares most here)

BulkFormer is the closest of the four models to toxicogenomics by modality: it is bulk, dense, magnitude-preserving, whole-genome, and it already fine-tunes on LINCS compound-perturbation and GDSC drug-response data — it has literally touched L1000-style perturbation transcriptomics. That is far nearer than scGPT (sparse single cells) or even Tahoe (single-cell cancer lines). But the residual gaps are exactly the toxicogenomics axes: pretraining is observational (no learned intervention prior), inputs are absolute log-TPM not fold-change (whereas DrugMatrix/TG-GATEs/L1000 speak in treated-vs-control log-fold-change), there is no dose or time conditioning, and it is human only (DrugMatrix/TG-GATEs are largely rat). So BulkFormer is the best starting point for a toxicogenomics FM and simultaneously a precise inventory of what such a model still needs — which §24.8 turns into a plan.

24.6 What is technically distinctive

Separating real advances from cheaper explanations (18 §18.7):

1. Dense-bulk modality. The core positioning: train on ≥14,000-gene dense samples, not sparse cells. A genuine niche, and the reason it beats single-cell FMs on bulk tasks — but that is partly a fair-matchup effect (single-cell FMs are out of their domain on bulk), so read the head-to-heads with that in mind.
2. Continuous expression + REE (magnitude kept). The most important representational choice, and a real improvement over binning/rank encoding for a modality where magnitude matters. This is scGPT §22.8-1 done natively.
3. Whole-transcriptome context via Performer. Keeping all 20,010 genes (no HVG subsetting) is only possible with linear attention — a real architectural enabler.
4. GCN co-expression prior fused per block. A hybrid inductive bias distinct from pure-attention FMs; helps, but also risks baking in correlation/confounds (§24.4, §24.7).
5. ESM2 protein-sequence gene initialization. A biological prior on gene identity, and (unused here) a latent route to cross-species transfer.

Genuine advance vs. cheaper explanation. Much of BulkFormer’s edge is (a) matching the model to the modality (bulk) and (b) magnitude-preserving continuous inputs — i.e. representation and data-fit, not a novel attention formula. Several head-to-head margins over single-cell FMs are thin (subtype 0.833 vs 0.830), and the biggest wins (imputation) are against baselines that co-expression already handles well. The honest read: a well-designed, well-matched bulk model, not a paradigm shift — but for toxicogenomics its representational choices are exactly the right ones.

24.7 Limitations and research gaps

Representation. Inputs are absolute log-TPM, not fold-change — the native currency of toxicogenomics (11_SC_FOUNDATION_MODELS.md §3) is treated-vs-control log-fold-change, which BulkFormer does not represent. No metadata at all (drug, dose, time, tissue) enters the model — a deliberate omission the authors flag. Protein-coding only, so regulatory non-coding RNAs (some of them toxicologically relevant) are invisible.

Architecture. The co-expression graph is fixed from raw Pearson correlations, which fold in batch/technical covariation; the GCN can therefore lean on non-biological edges. Bulk itself mixes cell types — a tissue average — and BulkFormer does not deconvolve, so cell-composition shifts (common under toxicant exposure — e.g. immune infiltration) are confounded into the sample vector.

Objective. Masked imputation is largely co-expression smoothing (§24.4) and is observational, so it teaches association, not dose-driven causation. The high imputation numbers should be read against strong co-expression baselines.

Evaluation. Several margins over single-cell FMs are thin, and single-cell FMs are out of domain on bulk (a favorable matchup). Prognosis is modest (AUROC 0.747). Train/val/test split details are under-documented, and there is no independent third-party replication of the core benchmarks yet.

Mechanism vs statistics. On the evidence, statistics: BulkFormer is an excellent model of co-expression structure in human bulk tissue. There is no demonstrated grasp of causal, dose-dependent perturbation mechanism — the perturbation task (LINCS PCC 0.493) is exactly where it is weakest.

Where it fails out-of-distribution. Rat / non-human (human vocabulary + human ESM2 init + human-only training — though the ESM2 route makes this the least hard of the failures to fix). Fold-change / treated-vs-control inputs (never seen). Dose and time (no representation). Interventional / toxicant exposure regimes at cytotoxic endpoints (observational pretraining). Non-coding transcripts. Cell-composition changes (no deconvolution). For toxicogenomics, the modality is right but every experimental axis is unmodeled.

24.8 How the model could be improved

The most detailed section, and unusually actionable here because BulkFormer is already the right modality — the improvements are mostly about adding the toxicogenomics axes, not rebuilding. Each: gap → how-to → data → metric+baseline → why it’s a contribution.

1. A fold-change (treated-vs-control) input/output mode

• Gap. Absolute log-TPM, but toxicogenomics is treated-vs-control log-fold-change (§24.7).
• How to. Add an input mode that ingests a paired (treated, control) profile and a head that outputs signed log-fold-change; reuse the REE continuous encoding for signed values.
• Data & metric. DrugMatrix / Open TG-GATEs / LINCS L1000 (all report fold-change vs control); evaluate signed-direction accuracy and magnitude correlation on held-out compounds against a per-gene mean-shift and a linear baseline.
• Why. It aligns the model’s currency with tox and makes it drop-in comparable to DrugMatrix/TG-GATEs analyses — the cousin of 22 §22.8-3 / 23 §23.8-5, but native to bulk.

2. Dose and time conditioning

• Gap. No dose/time axis (§24.7).
• How to. Add condition tokens (log-dose scalar + exposure-time token) to the currently metadata-free input; train with a monotonicity constraint in dose.
• Data & metric. Open TG-GATEs (dose × time in rat + human hepatocytes); hold out a dose/time and predict it; score dose-ordering and interpolation against linear-in-log-dose.
• Why. Concentration × time is the tox readout; this is the single most valuable addition, and TG-GATEs is purpose-built to evaluate it.

3. Cross-species (rat) via the ESM2 route it already has

• Gap. Human only — but DrugMatrix/TG-GATEs are largely rat (§24.7).
• How to. This is unusually cheap here: BulkFormer already initializes gene identity from ESM2 protein embeddings, so run rat genes’ proteins through the same ESM2 to get an ortholog-aligned rat gene vocabulary, then continue-pretrain on rat bulk. Orthologs land near their human counterparts by construction.
• Data & metric. Rat DrugMatrix/TG-GATEs; evaluate rat→human and human→rat transfer, stratified by ortholog conservation (20_REASONING_DATASET_IDEA.md §20.6).
• Why. It realizes the cross-species idea (22 §22.8-4) on a model that already has the mechanism half-built — the most tractable species bridge among the four models.

4. Interventional pretraining

• Gap. Observational GEO/ARCHS4 pretraining teaches association, not response.
• How to. Continue-pretrain (or co-train) on bulk perturbation corpora — LINCS L1000 (~1M+ bulk-like signatures), DrugMatrix, TG-GATEs — with the fold-change mode (#1) and dose/time tokens (#2), so the objective rewards predicting response.
• Data & metric. Held-out compounds/doses on L1000/TG-GATEs; the linear baseline (18 §18.4) as the honest bar.
• Why. It is the direct cure for the perturbation weakness (LINCS PCC 0.493) and turns a tissue-atlas model into a response model.

5. Metadata conditioning (compound, tissue, disease)

• Gap. Model ignores all metadata by design.
• How to. Re-introduce compound (SMILES), tissue, and disease tokens as conditioning, so the model can use context it currently discards; SMILES enables unseen-compound generalization (as in Tahoe-x1, 23).
• Data & metric. Tox and drug-response tasks with held-out compounds; measure unseen-compound response against a nearest-training-drug baseline.
• Why. Toxicogenomics is a metadata-conditioned prediction (compound, dose, time, tissue → response); a metadata-free model cannot do it.

6. Deconvolution-aware modeling

• Gap. Bulk mixes cell types; toxicant effects often are cell-composition shifts (immune infiltration, cell death) that a tissue average confounds.
• How to. Add an auxiliary cell-composition head (predict cell-type fractions, supervised by deconvolution references) so the model separates “the same cells changed” from “the mix of cells changed.”
• Data & metric. Tissues with paired bulk + single-cell references; evaluate composition recovery and whether conditioning on it improves response prediction.
• Why. It addresses a bulk-specific confound that directly distorts tox interpretation.

Cost note. #1, #2, #5, #6 are architecture/objective additions trainable by continuing from the released 147M checkpoint (single- to few-GPU with adapters); #3 reuses the existing ESM2 machinery; #4 is the heaviest (a perturbation pretraining run) but uses existing public corpora. Every one names its metric and baseline first (21 §21.11). Taken together, they are close to a blueprint for the toxicogenomics bulk foundation model the book keeps circling — which is precisely why BulkFormer, not scGPT or Tahoe, is the natural chassis to start from.

This chapter is exposition, not a result. BulkFormer facts are from Kang, Fan, Yi, Cui & Cui, “A large-scale foundation model for bulk transcriptomes,” bioRxiv 2025.06.11.659222 (Cell Systems 2026), with code on GitHub (KangBoming/BulkFormer) and weights on Zenodo (record 15559368): 522,769 human bulk samples from GEO+ARCHS4, 20,010 protein-coding genes, log(TPM+1) continuous inputs, ESM2 gene init + Rotary Expression Embedding, a GCN(co-expression) + Performer hybrid, ~147M main model (640-dim, 12-layer, 8-head), masked-value MSE at ~15%. Corrections to common assumptions flagged in-text: pretraining is GEO+ARCHS4 (not GTEx/TCGA, which are downstream), and the model is observational, human, protein-coding, absolute-expression, metadata-free. Uncertain items (⚠️): exact train/val/test split, the enumerated nine physiological systems, final Cell Systems volume/page, and the “150M”→147M-variant mapping. Reported metrics are developer-run with no independent replication yet, and several margins over single-cell FMs are thin. Next: Chapter 25, EVA — where we deliberately set the multimodal parts aside and ask what its transcriptomic representation, cross-species/-platform transfer, and scaling do that earlier single-cell FMs did not.