BulkFormer

Transcriptome foundation models have converged on single-cell RNA sequencing. Geneformer, scGPT, scFoundation, GeneCompass, and scLong are all pretrained on scRNA-seq, a modality that typically detects around 3,000 genes per cell. Bulk RNA-seq is a different measurement: roughly 16,000 genes per sample, dense rather than sparse, and the format in which most clinical cohorts, tissue atlases, and pharmacogenomic screens are actually recorded. Applying a single-cell model to bulk data means running it out of distribution, and it shows.

BulkFormer, from Qinghua Cui's group at Peking University with collaborators at Jilin University and Wuhan Sports University, is a foundation model built for the bulk modality from the start. It carries 150 million parameters over about 20,000 protein-coding genes and is pretrained on 522,769 human bulk transcriptomic profiles assembled from GEO and ARCHS4. Its encoder is a hybrid: a graph convolutional layer that propagates information along a prior gene–gene coexpression graph, followed by Performer layers whose linear-complexity attention makes a 20,000-token input tractable.

The preprint was posted in June 2025 and the peer-reviewed version appeared in Cell Systems in 2026. A December 2025 release withdrew the original compound-perturbation benchmark after the authors judged the LINCS assay platform too dissimilar from bulk RNA-seq to be a fair evaluation — an unusually candid retraction of one of their own headline results.

#Key Features

• Purpose-built for bulk RNA-seq: Pretrained on dense, whole-transcriptome samples rather than sparse single cells, closing a gap that left clinical and tissue-level tasks to models trained on the wrong modality.
• Hybrid GCN–Performer encoder: A graph convolution captures explicit gene–gene relationships from a coexpression prior while Performer attention models implicit global dependencies at linear cost.
• Protein-informed gene tokens: Each gene's initial representation is an ESM2 embedding of its canonical protein product, grounding gene identity in molecular function rather than an arbitrary index.
• Rotary expression embedding: Expression values are encoded as rotary positional representations, preserving both the ordering and the relative magnitude of expression levels without extra trainable parameters.
• One checkpoint, many tasks: The same frozen embeddings feed imputation, disease and cancer-subtype classification, prognosis, drug response, and gene essentiality prediction.

#Technical Details

Pretraining follows a BERT-style masked objective: 15% of gene expression values are masked and reconstructed under an MSE loss. Three input representations — ESM2-derived gene tokens, rotary expression embeddings, and an MLP-compressed sample-level vector — are summed. The coexpression graph retains the 20 highest-weight edges per gene and discards edges with Pearson correlation below 0.4. Downstream tasks take embeddings from the final Performer layer, max-pooled to a sample-level vector where required. Five checkpoints are released, from 37M to 147M parameters.

On held-out imputation with 15% masking, BulkFormer reached Pearson r = 0.954, against 0.806 for a variational autoencoder and below 0.150 for scFoundation and scLong, whose single-cell pretraining transfers poorly. Disease-type classification across 23 diseases from DiSignAtlas gave a weighted F1 of 0.939, ahead of scGPT at 0.885; on 33 TCGA cancer types it reached 0.833 against scGPT's 0.830. Survival prediction across roughly 10,000 TCGA patients gave AUROC 0.747 and AUPRC 0.549, against 0.726 and 0.520 for scFoundation. Drug sensitivity modeling on GDSC — 255 compounds across 700 cell lines — reached PCC 0.910 and SCC 0.879, and gene essentiality prediction on DepMap (17,862 genes, 1,103 cell lines) reached PCC 0.931 and SCC 0.759. All comparisons used ten-fold cross-validation.

#Applications

BulkFormer is aimed at the tissue- and patient-level analyses where bulk RNA-seq is the standard assay. Imputation recovers genes missing from older or shallowly sequenced datasets, which in the paper surfaced prognostic biomarkers previously obscured by dropout, including one whose high expression associated with a hazard ratio of 0.26. Its embeddings serve as features for disease and cancer-subtype classifiers, survival models, drug-sensitivity regression across cell line panels, and gene essentiality scoring that prioritizes cancer-specific vulnerabilities from expression alone. For groups working with TCGA, GTEx, or in-house clinical cohorts, it offers a pretrained representation where the alternative has been raw counts or PCA.

#Impact

BulkFormer fills a conspicuous hole in the transcriptome foundation model landscape by treating bulk RNA-seq as its own modality rather than a degenerate case of single-cell data, and it does so at a training cost the authors describe as substantially lower than that of the single-cell models it outperforms. Publication in Cell Systems and MIT-licensed weights on GitHub and Zenodo make it directly usable. The honest caveats are the authors' own: prognosis modeling remains difficult given the noise in bulk profiles, with an AUROC in the mid-0.7s leaving considerable headroom, and the withdrawn compound-perturbation benchmark is a reminder that a bulk-pretrained model does not transfer automatically to other expression assay platforms.