Illuminating cell states by a comprehensive and interpretable single cell foundation model

Domain: biology
Published: 2026-03-16
Editorial outcome: major_revision

Abstract
Advances in single-cell sequencing have enabled AI-driven foundation models with powerful data representation. However, their practical use is limited by real-world data sparsity, heterogeneity, and poor interpretability. To overcome these, we introduce ChrisCell. To enhance generalizability, we incorporate a large-scale single-cell dataset comprising 68 million cells, model parameters totaling 500 million, and challenging pretraining tasks. Notably, we introduce a Single-Cell Discretization (SCID) module that effectively represents cell embeddings, addressing data heterogeneity. For improved interpretability, the SCD module transforms high-dimensional and sparse single-cell data into a ”cell code,” facilitating recognition and analysis. Additionally, we also present ChrisCell-Graph, a plug-and-play tool that integrates ChrisCell’s features with multimodal data (genes, cell communication, annotations) to build a knowledge graph for biological discovery. Extensively evaluated, ChrisCell outperforms strong baselines in all downstream tasks, and also uncovered intriguing biological phenomena with compelling explanations. ChrisCell aspires to serve as a truly applicable and generalizable AI tool for the cell biology community.

Introduction
Cells, as the fundamental units of life, though tiny, execute a remarkable diversity of essential functions crucial for maintaining human health [1, 2]. The discovery of cells by Robert Hooke through his microscope captivated scientists, prompting an enduring fascination with their complexity and intricacies, and initiating a quest to uncover the secrets of cellular life. Historically, during the microscopy era, biologists primarily speculated about cells at the microimaging level, where observations were largely confined to their morphology, activity, and interactions. This limitation impeded deeper investigations into the complexities of cellular states. However, with the advent of the post-genomic era, a paradigm shift has emerged toward a molecular understanding of cell states, highlighting the myriad factors that influence them [3]. A significant international initiative, the Human Cell Atlas Project [4], was launched in 2017 to establish ”an open, comprehensive reference map of the molecular state of cells in healthy human tissues.” This project employs single-cell RNA sequencing (scRNA-seq)[5, 6] techniques and has inspired similar endeavors focusing on well-studied animal models[7, 8] and patient-derived tissue samples [9] across a diverse range of diseases.

In addition to RNA sequencing, other omics technologies, such as single-cell ATAC sequencing (scATAC-seq) [10, 11] for chromatin accessibility and proteomics, have been employed to comprehensively analyze cell states. Concurrently, the availability of sequenced single-cell data has surged to tens of millions, presenting a significant opportunity for the application of AI-based methods in cell biology.

Leveraging large-scale scRNA-seq data, several deep learning-based methods have emerged to learn informative cell representations through pretraining on extensive datasets. For instance, scBNERT[12] employs a bidirectional transformer architecture[13] to pretrain on unlabelled scRNA-seq data, subsequently fine-tuning for specific tasks. Similarly, scGPT[14 ] utilizes a generative pretrained transformer trained on a repository of over 33 million cells, demonstrating robust performance in both generative and representation tasks. scFoundation[15] adopts an asymmetric encoder-decoder transformer to conduct read-depth-aware modeling pretraining on 50 million transcriptomic profiles. Additionally, scPRINT [16] focuses on gene network inference by interpreting the attention matrix generated by its transformer, which is pretrained on 50 million cells. Collectively, these models leverage large-scale single-cell data to learn robust cell representations and exhibit impressive performance across various downstream tasks. Despite these advancements, emerging evidence indicates notable generalizability limitations of these foundation models [17, 18]. These limitations can be categorized into three primary aspects: (1) Sparsity of Single-Cell Data [19]: In scRNA-seq and scATAC-seq datasets, only a subset of genes or regions is expressed or captured, resulting in a substantial amount of zero gene expression. This sparsity challenges the model’s capacity to extract compact and rich information from high-dimensional data. (2) Heterogeneity of Single-Cell Data [17, 18]: Single-cell data are generated by various institutes using different sequencing, processing, and normalization methods. This diversity leads to significant variations in the sequenced genes and expression scales, complicating the generalization of single-cell foundation models. (3) Poor interpretability of Foundation Models: These models contain a vast number of parameters and exhibit high complexity. While they are capable of generating informative gene and cell representations, the high-dimensional and continuous nature of these representations renders them challenging for biological scientists to comprehend.

We present ChrisCell, a comprehensive and interpretable single-cell foundation model. ChrisCell innovatively integrates a Single-Cell Discretization (SCD) module into the single-cell foundation model. This module utilizes a unified cell codebook to transform the cell representation into a cell code and uses an SCD cell embedding derived from the cell codebook to represent the cell corresponding to the code. ChrisCell employs an encoder-SCD-decoder architecture, encompassing 511 million parameters and pretrained on over 68 million single-cell data points. The model incorporates multiple pretraining tasks at both the cell and gene levels. At the cell level, the task involves predicting cell properties such as type, tissue, age, sex, and disease. At the gene level, we implement a cell-gene-prompt (CGP)-based Zero-Inflated Negative Binomial (ZINB) prediction task [20, 21], which effectively enhances cell representation and captures gene co-expression patterns, and addresses the high sparsity of the data. Notably, our findings and prior evidence [22, 23] indicate that the SCD module effectively mitigates data heterogeneity and batch effects [24]. By utilizing the unified representation and cell codes from the cell codebook, the SCD module represents cells from different batches, leading to significant enhancements in both generalization and interpretability. Additionally, we introduce a lightweight model, ChrisCell-Graph, designed to integrate multi-omic and multimodal data to construct a Graph Knowledge Network (GKN). By leveraging the representations generated by ChrisCell, the plug-and-play functionality of ChrisCell-Graph facilitates the transfer of insights to diverse real-world datasets, enabling the discovery of biological knowledge through the GKN.

We evaluated the performance of our models across multiple tasks, categorized into prediction and analysis tasks. In the realm of prediction tasks, which encompass cell clustering, annotation, property prediction, gene and drug perturbation prediction, as well as drug response prediction, ChrisCell demonstrates significant advancements over other foundation models such as scGPT and scFoundation. For analysis tasks, we integrated ChrisCell and ChrisCell-Graph into four distinct tasks: cell state discovery, gene discovery, gene regulatory network (GRN) analysis [25, 26], and multimodal analysis. These analyses highlight the unique advantages of our models in terms of interpretability. We contend that the introduction of ChrisCell and ChrisCell-Graph effectively bridges the gap between AI models and biological research, thereby enhancing the understanding of cellular and genetic mechanisms.

Results
2.1 Overview of ChrisCell
Fig. 1 illustrates the comprehensive framework of ChrisCell. Fig. 1a depicts the model architecture of ChrisCell. For the encoder input, we first append a [Count] token, representing the sum of gene expression values. We then extract expression values from gene expression profiles and randomly sample genes as input for the encoder. The sampling rates follow a uniform distribution U(0.5, 0.8) to facilitate multi-scale information extraction. These sampled gene expressions are processed through an embedding module comprising three components: gene expression embedding, gene position (gene id) embedding, and gene ontology embedding. The gene expression embedding transforms continuous gene expression scalars into learnable high-dimensional vectors, while the gene position embedding encodes gene identification(gene names, like DLK1, CLU ). The Gene Ontology embedding module incorporates prior knowledge from Gene Ontology [27] to model gene co-expression, which is updated during training. Subsequently, these embeddings are passed through 12 Transformer layers to refine the gene embeddings, leading to the generation of a cell embedding via a pooling method. We input the gene embeddings into an SCD module to obtain SCD gene embeddings, while the cell embedding is quantized to produce an SCD cell embedding and a cell code, which serves as a data-driven identifier for the cell. The SCD cell embedding is utilized for multiple property predictions to train both the encoder and the SCD module. In the decoder component, we initialize all expression values to zero, allowing the embedding module to focus solely on capturing gene position and ontology information. We replace the decoder embeddings at the corresponding gene positions with SCD gene embeddings as prompts, setting the SCD cell embedding at the ”Count” token for the cell prompt. These processed decoder embeddings are subsequently input into six Performer [28] layers and a Zero-Inflated Negative Binomial (ZINB) predictor to estimate the ZINB distribution for each gene expression value. We select predictions for non-sampled genes and a subset of non-expressed genes (equal to or fewer than the non-sampled) to compute the ZINB loss against the corresponding gene expression values, facilitating the training of the entire model. Further architectural details are provided in Section 4.

Fig. 1b visualizes the six prediction tasks implemented using ChrisCell, including cell clustering, annotation, property prediction, gene perturbation prediction, and drug response prediction. All prediction tasks were performed under the zero-shot setting without fine-tuning. Additionally, we provide the comprehensive ablation study of ChrisCell in Supplementary Section 1; and the results of batch correction tasks to demonstrate that ChrisCell effectively resolves data heterogeneity in Supplementary Section 3. We also showcase the extensive experiments of ChrisCell on Cell clustering and cell type annotation tasks, comparing with broader baselines on TS2 [29, 30] datasets in Supplementary Section 2.

Fig. 1c illustrates the integration of ChrisCell and ChrisCell-Graph. ChrisCell-Graph receives representations from ChrisCell and combines these with other omics and modal data, such as scATAC-seq and property annotation data, to construct the graph knowledge network. Fig. 1d presents four analysis tasks: cell state discovery, gene discovery, GRN analysis, and multimodal analysis. In all analysis tasks, we only utilize the embeddings provided by ChrisCell pre-training model and train ChrisCell-Graph on the corresponding dataset. Comprehensive details regarding the datasets and experiments are available in the Supplementary Sections 12 and 13.

2.2 Representation improvement on cell clustering
The cell clustering task is fundamental to single-cell analysis. By clustering cells based on similar gene expression profiles, biologists can achieve cell type identification, analyze cell heterogeneity, conduct drug response studies, and discover biomarkers. Additionally, the effectiveness of the clustering task serves as a rigorous test of representation quality; superior cell representations yield higher accuracy and clearer visualizations.

To evaluate this, we applied ChrisCell to three scRNA-seq datasets, specifically the Pancreas [14], Zheng 68K [31], and Segerstolpe [32] datasets, all of which include cell-type annotation labels. We compared results with raw gene expression, GeneCompass [33], scBERT [12], UCE [34], CellFM [35], scGPT, and scFoundation. For each method, we derived cell representations and employed the K-means algorithm for clustering based on cell type labels. We utilized several evaluation metrics: Accuracy, Normalized Mutual Information [36] (NMI), Adjusted Rand Index [37] (ARI), and Silhouette Coefficient [38] (SIL). Among these metrics, Accuracy, NMI, and ARI assess whether cells of the same type are clustered together, while SIL measures the closeness of samples to their respective clusters and their separation from other clusters. Thus, these metrics provide a comprehensive evaluation of ChrisCell’s performance in clustering tasks.

Fig. 2a displays the overall performance of the eight methods across the three datasets. The metrics for ChrisCell on the Pancreas dataset are as follows: NMI (0.9424), ARI (0.9492), and SIL (0.4636). For the Zheng dataset, the metrics are: Accuracy (0.5127), NMI (0.6000), ARI (0.3904), and SIL (0.2050). The Segerstolpe dataset yields: Accuracy (0.7401), NMI (0.8720), ARI (0.6824), and SIL (0.6288). ChrisCell achieves exceptional performance across all datasets and metrics, particularly in the Pancreas dataset, where Accuracy, NMI, and ARI exceed 0.94, with a SIL of 0.4636. In contrast, other methods report NMI below 0.89, ARI below 0.77, and SIL below 0.34. Furthermore, in the Segerstolpe dataset, ChrisCell outperforms other methods across all metrics, leading by 0.23 in NMI, 0.30 in ARI, and 0.5 in SIL. These results underscore the superior representation and zero-shot capabilities of ChrisCell.

Fig. 2b visualizes UMAP [39] representations with cell type labels for the Segerstolpe dataset based on the eight methods. The UMAP representation of raw gene expression appears scattered, indicating that the high dimensionality and sparsity of raw gene expression inhibit the capture of critical information. In contrast, other methods, like scGPT and scFoundation achieve more precise clustering of cells of the same type; however, they still fail to assign all cells into correct and compact clusters. Notably, ChrisCell successfully separates all cell types and compresses each cluster into a compact group, effectively capturing the biological relationships between different cell types. For instance, endothelial cells and pluripotent stem cells (PSCs) are positioned closely in the UMAP generated by ChrisCell. This proximity reflects the biological understanding that PSCs can differentiate into endothelial cells and co-express certain marker genes, such as VEGF (Vascular Endothelial Growth Factor), CD31, and CD34, which are commonly used to label endothelial and precursor cells. Consequently, ChrisCell adeptly identifies the underlying relationship between these two cell types while exhibiting similar representations.

2.3 Cell annotation and properties prediction
Cell type annotation is a critical step in the analysis of single-cell data, aimed at assigning biological significance to clustered cell populations by correlating them with known biological cell types or states, as defined by established marker genes. To further validate the zero-shot capability and representation quality of ChrisCell, we applied it to three scRNA-seq datasets, similar to those used in the clustering tasks, and compared the results with four baseline methods. We established a linear classification model using scikit-learn [40], employing the representations generated by the eight methods as input to predict cell types. The evaluation of zero-shot performance utilized Accuracy, F1 score, Precision, and Recall as metrics.

Fig. 2c presents the overall performance of the eight methods across the three datasets. The metrics for ChrisCell are as follows: Pancreas (Accuracy: 0.88, F1 score: 0.87, Precision: 0.88, Recall: 0.88), Zheng (Accuracy: 0.90, F1 score: 0.89, Precision: 0.90, Recall: 0.89), and Segerstolpe (Accuracy: 0.85, F1 score: 0.80, Precision: 0.76, Recall: 0.85). Consistent with the clustering tasks, ChrisCell demonstrates significant superiority across all metrics for the three datasets. Specifically, ChrisCell outperforms the other methods by margins of 0.20, 0.22, 0.21, and 0.22 in the Pancreas dataset, 0.26, 0.27, 0.26, and 0.25 in the Zheng dataset, and 0.26, 0.33, 0.26, and 0.26 in the Segerstolpe dataset across the four metrics. These results indicate that the representations provided by ChrisCell are more informative and exhibit superior linear interpretability.

To evaluate ChrisCell’s performance on multiple property prediction tasks, we created a test dataset that excluded samples from the pretraining dataset and included cell property annotation labels such as age, disease, sex, and tissue type. We applied the eight methods to this dataset and assessed performance using accuracy as the metric.

Fig. 2d showcases the accuracy of the eight methods on multiple property prediction tasks. Notably, the accuracy of scGPT and scFoundation on these tasks is lower than that of raw gene expression, indicating their inability to capture critical patterns related to these properties, resulting in information loss. In contrast, ChrisCell achieves higher accuracy in property predictions, demonstrating its capacity to extract comprehensive insights and transform raw gene expression values into a compact yet informative cell representation.

Fig. 2e illustrates the UMAP visualization of ChrisCell on this test dataset, with tissue type labels. The results indicate that ChrisCell effectively aggregates cells of the same tissue type into cohesive populations. While some groups of identical tissue types are separated, these distinctions reflect different cell types or states, underscoring that ChrisCell provides a more nuanced understanding of cellular patterns.

2.4 ChrisCell enhances bulk-level RNA-seq data representation
Bulk RNA-seq [41] data is generated by extracting RNA from a mixed population of cells, rather than from individual cells, and subsequently sequencing it. This approach provides the average gene expression profile of all cells in the sample, in contrast to single-cell data, which offers single-cell resolution. Although bulk RNA-seq data lacks the cellular heterogeneity captured by single-cell sequencing, it remains essential for numerous practical applications, such as tumor and drug response research, due to its more established and cost-effective experimental and analytical processes. To evaluate whether ChrisCell can enhance the representation of bulk RNA-seq data and further demonstrate its generalization capabilities, we applied ChrisCell to two drug response prediction tasks: drug response classification and cancer drug response prediction.

The goal of drug response classification is to predict the sensitivity of a patient or cell to a specific drug based on various biological characteristics, such as gene expression. This task offers valuable insights for precision medicine, optimizing drug development, and identifying potential new therapies. Our approach incorporates both bulk and single-cell data as inputs and outputs the predicted sensitivity for each cell. Utilizing the downstream module SCAD [42], we leverages ChrisCell to obtain the embeddings of single-cell and bulked data to train SCAD models for predicting a cell’s sensitivity to the drug. We compare ChrisCell performance with that of two foundation models: scGPT and scFoundation under the same setting. We also include Linear and Random Forest [43] as additional baselines for more comprehensive comparison. Our analysis focused on four drug test sets which are collected from the Genomics of Cancer Drug Sensitivity database [44]: Sorafenib, PLX4720-451Lu, NVP-TAE684, and Etoposide. We measured the quality of predictions using the AUC and recall metrics.

Fig. 3a presents the ROC curves for the six methods across four drugs. ChrisCell demonstrates consistently superior performance, particularly for Etoposide, where it achieves higher accuracy than other baselines at all thresholds. Notably, the baseline SCAD, Linear and Random Forest exhibit poor results for PLX4720 and Sorafenib, performing worse than random predictions, indicating significant challenges in predicting sensitivity for these drugs. In contrast, ChrisCell maintains a clear advantage, especially with Sorafenib, where the true positive rate exceeds 0.8 when the threshold is above 0.2, highlighting its effectiveness.

Fig. 3b illustrates the average AUC and recall results for the six methods across the four drugs. The average AUC and recall for ChrisCell are as follows: Sorafenib (0.88, 0.91), PLX4720-451Lu (0.77, 0.80), NVP-TAE684 (0.83, 0.86), and Etoposide (0.76, 0.79). These results indicate that ChrisCell maintains both AUC and recall above 0.75, achieving leading performance across all drugs and metrics compared to other methods. This demonstrates that ChrisCell effectively transfers to bulk RNA-seq data, extracting more information to enhance the accuracy of drug response classification.

In addition to drug response classification, we applied ChrisCell to the cancer drug response (CDR) prediction task, which examines tumor cell responses to drug interventions. Given the drug information and bulk RNA-seq data, the model predicts the half-maximal inhibitory concentration IC50. To enhance prediction quality, we utilized ChrisCell to encode the bulk data, obtaining cell representations, and subsequently employed the downstream module DeepCDR [45] for the CDR prediction task. The prediction quality was assessed using Pearson correlation coefficients.

Fig. 3c compares the results of ChrisCell with other baselines across all cell lines, cancer types, and drug data. These results demonstrate that ChrisCell achieves stronger Pearson correlation coefficients than other baselines across the three test sets. Furthermore, Fig. 3d indicates that ChrisCell attains the highest average Pearson correlation coefficients of 0.923, compared to 0.780 for DeepCDR, 0.852 for scGPT, and 0.878 for scFoundation. Fig. 3e displays the overall distribution of Pearson correlation coefficients for the four methods. ChrisCell not only achieves the highest average Pearson correlation coefficients but also maintains stable variance across individual samples. These results provide further evidence that ChrisCell possesses significant capability in transferring to bulk RNA-seq data and generalizing to unseen datasets and tasks.

2.5 ChrisCell enhance gene representation for perturb-seq data
Perturb-seq [46] is a high-throughput experimental technique that integrates CRISPR-based gene perturbation [47] with scRNA-seq. This method enables researchers to systematically investigate the effects of gene perturbations on gene expression at the single-cell level, thereby facilitating a deeper understanding of cellular functions, gene interactions, and biological pathways.

To evaluate whether ChrisCell enhances gene-level representation and captures gene co-expression during pretraining, we applied ChrisCell to three perturb-seq datasets: Adamson [48], Norman [49], and Dixit [46], focusing on gene perturbation prediction tasks. We combined ChrisCell with an advanced model known as GEARS [50] for predicting single-cell-resolution perturbations. The original GEARS model utilized a Gene Ontology knowledge graph to represent unseen gene perturbations by learning from a combination of previously observed gene perturbation nodes. Additionally, it employed a gene co-expression graph, integrated with perturbation information, to predict post-perturbation gene expression.

We compared ChrisCell with GEARS, Linear, GeneCompass, CellFM, scGPT, and scFoundation, all under the same experimental conditions. The evaluation of prediction accuracy was conducted using several metrics: Mean Squared Error (MSE), MSE DE (representing the distance between the top 20 differential gene expressions of the ground truth and corresponding predictions), Pearson correlation coefficients, Pearson correlation coefficients DE, and Pearson delta, which calculates the Pearson correlation between the expression changes of prediction and ground truth. MSE measures the overall distance between predictions and ground truth, while Pearson correlation coefficients and Pearson delta assess the correlation between predicted and observed post-perturbation expression changes.

Fig. 4a-d presents the overall performance of ChrisCell and other baseline models across three perturb-seq datasets, evaluated using four metrics: MSE, MSE DE, Pearson correlation coefficients, and Recall. The performance metrics for ChrisCell on the three datasets are as follows: Adamson (0.032, 0.830, 0.990, 0.850), Norman (0.020, 0.200, 0.993, 0.900), and Dixit (0.013, 0.082, 0.990, 0.970). These results indicate that ChrisCell outperforms all other models across all datasets and metrics, demonstrating a clear advantage. Notably, the MSE and MSE DE values for ChrisCell on the three datasets are consistently below 0.05 and 0.85, respectively, indicating that ChrisCell’s predictions are substantially close to the ground truth for both overall genes and the top 20 differentially expressed genes. Furthermore, the Pearson correlation coefficients (PCC) and PCC DE values for ChrisCell exceed 0.99 and 0.85, respectively, underscoring the substantial similarity of its predictions to the ground truth. In addition, we observe that the linear model performs comparably to some deep learning models. Notably, on the Dixit dataset, the MSE, MSE DE, PCC, and PCC DE are 0.063, 0.174, 0.82, and 0.74, respectively. These results are comparable to, or even better than, those of GEARS (0.145, 0.315, 0.68, 0.53), GeneCompass (0.084, 0.182, 0.7, 0.63), CellFM (0.07, 0.14, 0.85, 0.74), scGPT (0.076, 0.11, 0.73, 0.70), and scFoundation (0.066, 0.15, 0.83, 0.74). However, on the Adamson and Norman datasets, the linear model achieves performance that is only comparable to or lower than that of GEARS, where these two datasets contain more gene perturbation data. This indicates that while the linear model performs well on smaller datasets, foundation models tend to generalize better on larger-scale datasets.

Fig. 4e illustrates the Pearson delta comparison of ChrisCell with other baselines across the three datasets. The figures demonstrate that ChrisCell predicts a greater number of samples with higher Pearson delta values than the other baselines. The average Pearson deltas for ChrisCell are 0.631, 0.451, and 0.444 for the Adamson, Norman, and Dixit datasets, respectively. In comparison, the Pearson deltas for the other baselines are as follows: GEARS (0.050, 0.189, 0.090), CellFM(0.483, 0.376, 0.243) scGPT (0.464, 0.258, 0.236), and scFoundation (0.564, 0.388, 0.256). These results indicate that ChrisCell achieves superior prediction performance, not only for individual samples but also in terms of overall Pearson delta, highlighting the effectiveness of ChrisCell’s gene representation and its capability to capture gene co-expression.

Fig. 4f showcases the expression changes of the top 20 genes following the two-gene perturbation with LYL1 and CEBPB. ChrisCell’s predictions exhibit closer and more consistent delta values compared to the control for almost all genes, indicating that the individual gene representations produced by ChrisCell are robust and informative.

2.6 ChrisCell enhance gene representation for drug perturbation prediction
Drug perturbation prediction aims to elucidate the effects of specific drugs on cellular gene expression, metabolic activity, and other biological characteristics, with a wide range of practical applications that have propelled advances in medical and pharmaceutical research. By predicting the effects of drugs on cells, researchers can identify new drug targets and anticipate an individual’s response to specific medications based on genomic and biomarker profiles, thereby enhancing treatment efficacy.

To evaluate whether ChrisCell’s gene representation can improve drug perturbation prediction, we applied ChrisCell to a single-cell perturbation dataset derived from human peripheral blood mononuclear cells (PBMCs), referred to as OP3 [43]. Fig. 5a illustrates the pipeline for conducting drug perturbation prediction on the OP3 dataset. We utilized a Uni-Mol v1 pretrained model [51] to encode the drug SMILES data, thereby obtaining the drug representation. Concurrently, we employed ChrisCell to derive gene representations from the provided scRNA-seq data. These gene embeddings were then concatenated with the drug representation and input into a multi-layer perceptron to predict differential gene expression. We performed this prediction using various methods, including raw gene expression data, scGPT, and scFoundation, under consistent settings. In addition, to more comprehensive comparison, we also include non-foundation models, like scVIDR [52] and CellOT [53]. The evaluation metrics included MSE, MSE DE, Pearson correlation coefficients(PCC), Pearson correlation coefficients DE(PCC DE), and Pearson delta, with Pearson delta DE serving as an additional evaluation metric.

Fig. 5b presents the overall performance of ChrisCell and other baseline methods on the OP3 datasets, evaluated using four metrics. The results indicate that relying solely on raw gene expression data yields suboptimal performance, with MSE and MSE DE values exceeding 1.25 and 4, respectively, while the PCC and PCC DE values fall below 0.66 and 0.46. This suggests that predicting using raw gene expression data does not adequately align with the ground truth. Besides, we noticed that the non-foundation models can achieve comparable or better performance over some foundation models. For instance, the MSE, MSE DE, PCC and PCC DE of CellOT are 0.247, 0.763, 0.901 and 0.782, which are superior than scGPT(0.887, 2.333, 0.853, 0.678) and scFoundation(0.312, 1.163, 0.897, 0.800). However, ChrisCell still achieves better performance over other methods. The representation derived from ChrisCell reduces MSE and MSE DE to 0.087 and 0.233, respectively, while improving PCC and PCC DE to 0.965 and 0.945. These results demonstrate a significant enhancement in prediction quality, with ChrisCell outperforming other robust baselines, such as scVIDR and scFoundation, across all metrics.

Fig. 5c compares the Pearson delta and Pearson delta DE of ChrisCell with five baseline methods on the OP3 dataset. The average Pearson delta and Pearson delta DE for ChrisCell are 0.478 and 0.677, respectively. In comparison, raw gene expression yields values of 0.12 and 0.223, while scGPT reports 0.235 and 0.356, scVIDR is 0.212, 0.343, CellOT reports 0.268, 0.391, and scFoundation yields 0.267 and 0.388. ChrisCell exhibits superior performance in both Pearson delta metrics, indicating that its gene representation captures more co-expression information, thereby facilitating more accurate predictions.

Fig. 5d showcases the expression changes of the top 20 genes following drug perturbation with Clotrimazole in NK cells. The predictions from ChrisCell demonstrate closer and more consistent delta values compared to the control across nearly all genes, indicating that the individual gene representation produced by ChrisCell is both robust and informative.

2.7 Combining ChrisCell and ChrisCell-Graph on analysis tasks
Fig. 6a visualizes the pipeline for integrating ChrisCell and ChrisCell-Graph. Given the molecular census data of a single cell, ChrisCell generates both cell and gene embeddings for ChrisCell-Graph. ChrisCell-Graph constructs a knowledge graph in which each cell, molecule, and property is represented as a node within a unified framework, connecting all factors associated with the cell. This approach facilitates the embedding of these properties into a shared latent space. Subsequently, ChrisCell-Graph employs a graph attention module [54] to update cell representations and compute attention weights for each neighboring node. This module is designed to incorporate multi-modal information into the cell representation, thereby enhancing the interpretability of the model.

2.8 Marker gene discover and analysis with ChrisCell
Fig. 6b presents the UMAP visualization of raw gene expression and ChrisCell embeddings, labeled by cell types in the pancreas dataset. The UMAP representation of raw gene expression fails to distinctly separate the various cell types, particularly for smooth muscle, endothelial, and myeloid cells. In contrast, ChrisCell effectively segregates these cells into distinct populations based on their respective cell types, demonstrating the interpretability of the ChrisCell representation. We utilized ChrisCell to encode the scRNA-seq data from the pancreas dataset and employed max pooling to derive cell embeddings. An importance metric was defined to represent the significance of each gene for the corresponding cell types, calculated as the frequency with which a gene’s channel is extracted during max pooling. This frequency was normalized by dividing the number of times each gene is extracted by its maximum value.

Fig. 6c and Fig. 6d illustrate the importance of the top 20 genes in Acinar and Endocrine cells, as predicted by ChrisCell. Based on the importance rankings provided by ChrisCell, we successfully identified several marker genes for both Acinar and Endocrine cells. For instance, the top 10 genes, including CEL [55], CTRB2 [56], CPA2 [57], CRTB1 [56], SPINK1 [58], CLPS [59], and SYCN [60], are recognized as marker genes for Acinar cells. These genes are highly expressed in Acinar cells and are primarily involved in the synthesis and secretion of digestive enzymes, lipid metabolism, and cytoprotective mechanisms. Notably, CTRB2 [56] and CRTB1 [56] act as precursors for digestive enzymes, while CPA2 [57] and CEL [55] contribute to protein and fat digestion, working synergistically to ensure adequate availability of digestive enzymes in the gut. In the case of Endocrine cells, the top 20 genes identified by ChrisCell, including CHGA [61], INS [62], FTL [63], SST [64], and GCG [65], serve as crucial marker genes. These genes are integral to hormone synthesis, storage, and secretion, thereby influencing metabolic and physiological homeostasis throughout the body. Additionally, genes such as CADM1 [66] play significant roles in cell aggregation and interactions, contributing to the structural integrity and functionality of pancreatic endocrine cells by promoting intercellular adhesion, which affects hormone secretion and cell growth.

Fig. 6e displays the top 100 gene importance scores for pancreatic ductal cells. ChrisCell not only highlights the significance of important genes or marker genes, such as KRT19 [67], SOX9 [68], CDH1 [69], and PDX1 [70], but also captures unique markers for specific subtypes of ductal cells, including HMMR [71] and TPX2 [72]. Furthermore, ChrisCell emphasizes a wide range of genes whose functions are essential for myeloid cells.

2.9 Cell state discover and analysis with ChrisCell-Graph
To evaluate the efficacy of ChrisCell-Graph in analyzing scRNA-seq data for cell state discover, we conducted an experimental analysis using a well-established brain dataset. Fig. 7a and b illustrate the UMAP visualization of the brain dataset, depicting both raw gene expression and the ChrisCell-Graph representation, respectively. The ChrisCell-Graph representation distinctly separates major cell categories, such as neurons, astrocytes, central nervous system macrophages, oligodendrocytes, and oligodendrocyte precursor cells, offering clearer differentiation compared to the raw gene expression. Moreover, ChrisCell-Graph effectively disentangles smaller cell categories, including leukocytes, endothelial cells, ependymal cells, and astrocytes. We also provide further analytical results of ChrisCell-Graph on cell-type-specific datasets, along with ablation studies, in silico perturbation experiments, integration with other foundation models, and a comparison of GRN analysis with CollecTRI [73] in Supplementary Sections 4–8, respectively.

Fig. 7c presents the UMAP visualization of the ChrisCell-Graph representation, with labels corresponding to the cell states. This representation enables a more precise state identification of cells based on known cell type categories. For Fig. 7d, we obtain the corresponding SCD embeddings of each cell, performing the UMAP visualization. Fig. 7e displays a heatmap illustrating the relationship between cells with different cell states and genes, utilizing the average attention weights derived from the graph attention module.

According to the heatmap in Fig. 7e, we identify crucial genes with high attention weights for each cell state, which are annotated in Fig. 7d and annotated some state genes in Fig. 7c for the following analysis examples. We denote these genes as state genes. This figure showcases the UMAP visualization based on the SCD cell embeddings, annotated with these significant genes. Subsequent analyses reveal several notable findings related to cell states.

According to the results from Fig. 7c and d, we argue that ChrisCell-Graph’s characterization assists biological researchers in swiftly identifying crucial genes essential for understanding cell states. For instance, by comparing Fig. 7b and c, it exhibits that cells categorized under state 2289 predominantly fall within neuron and oligodendrocyte precursor cell (OPC) types, with RBFOX1 and ANKS1B identified as high-attention-weight genes. Notably, both genes are recognized as vital for neuronal function. RBFOX1 [74] (RNA binding protein, Fox-1 homolog 1) is an RNA-binding protein involved in splicing regulation, playing a critical role in neuronal development and function by regulating neuron-specific gene splicing [74]. ANKS1B [75] (ankyrin repeat and SAM domain-containing protein 1B) contributes to cell signaling and intercellular interactions, influencing synapse formation and plasticity, and is associated with various neurological disorders, including autism spectrum disorders, thereby impacting neurodevelopment and function. Furthermore, these genes also provide essential functions for oligodendrocyte precursor cells. In OPCs, RBFOX1 may influence RNA splicing and post-transcriptional regulation, affecting their proliferation and differentiation. ANKS1B is involved in cell signaling and neural development, particularly in interactions with neurons. Consequently, ChrisCell-Graph classifies neurons and OPCs under the same cell state, with corresponding proximity in the visualization of Fig. 7c.

The second example involves three distinct states of oligodendrocytes, identified by states 3430, 3302, and 3174 by comparing Fig. 7b and c. We also annotate state genes corresponding to these three states in Fig. 7c by referring to Fig 7d. In the cells corresponding to state 3430, the PCDH9 [76] gene exhibits a significant attention weight. PCDH9 [76] (Protocadherin 9) is recognized as a critical gene in oligodendrocytes, belonging to the protocadherin family of proteins that facilitate cell-cell adhesion and interactions, particularly within the nervous system. The expression of PCDH9 is closely associated with oligodendrocyte maturation and myelin formation, indicating that cells with state 3430 are likely in a mature oligodendrocyte state [77]. In the cells designated by state 3302, we observe high attention weights for the ANK3 [78] and FTX [79] genes. FTX [79] (Ftx RNA) is a long non-coding RNA primarily involved in the nervous system and oligodendrocytes; it may regulate myelin formation as well as cell proliferation and differentiation. ANK3 [78] (Ankyrin 3) is a critical cytoskeletal protein that plays a central role in the nervous system, contributing to the construction and maintenance of the cytoskeleton. It connects the cell membrane to intracellular structures and is expressed in both neurons and glial cells, facilitating signaling and intercellular interactions with a closer relationship to neurons. The prominence of these two genes in state 3302 suggests that these cells are engaged in myelin formation and neural development, closely associated with neuronal activity [80]. Notably, Fig. 7c shows that state 3302 cells are positioned closer to neurons and oligodendrocyte precursor cells (OPCs) in the UMAP visualization compared to the other two oligodendrocyte states. For cells classified under state 3174, both PCDH9 [76] and ANK3 [78] exhibit prominent expression. Given that ANK3 is a key gene in state 3302 and PCDH9 is crucial in state 3430, this finding suggests that cells in state 3174 represent a transitional state between states 3302 and 3430. This classification allows us to differentiate a broad group of oligodendrocytes into three distinct and meaningful states. Fig. 7f and g provide further validation of our analysis. In terms of the PCDH9 gene, cells in state 3430 show the highest attention weight, while cells in state 3174 also exhibit a relatively high attention weight, albeit lower than that of state 3430. Conversely, state 3302 cells show minimal attention weight for the PCDH9 gene. Regarding the FTX gene, most cells in state 3302 display relatively high attention weights, while state 3174 cells possess somewhat lower weights. In contrast, state 3430 cells do not exhibit significant attention weights for the FTX gene, with many not even considering its weight.

The third example involves two distinct states, 1490 and 1474, within central nervous system macrophages in Fig. 7b and c. In cells classified under state 1474, the DOCK4 [81] gene shows high attention weight, whereas in state 1490 cells, PAG1 [82] exhibits a similar prominence. Although both DOCK4 and PAG1 are highly expressed in immune cells, they fulfill different functions depending on the context. PAG1 [82] (Phosphoprotein Associated with Glycosphingolipid-enriched Microdomains 1) is a membrane protein that modulates cellular signaling by regulating tyrosine kinase activity, playing a crucial role in immune cell activation and signaling, particularly within T and B cells. In contrast, DOCK4 [81] (Dedicator of Cytokinesis 4) regulates small GTPases (e.g., Rac1) and orchestrates cytoskeletal reorganization and cell migration through the activation of these GTPases, thus playing a key role in cell motility, development, and immune responses. The functional divergence between these two genes indicates that the primary roles of the cells associated with the two codes differ, with state 1490 cells actively engaging in immune functions [82], while state 1474 cells are primarily involved in migration [81].

These findings and examples indicate that ChrisCell-Graph is capable of learning the intricate relationships and differences between cells in a data-driven manner, effectively identifying discrete cell states using the SCD module. Moreover, the incorporation of the attention module effectively enhances the interpretability of each cell state, highlighting crucial genes that are pivotal for recognizing the characterization of specific cell states. By providing a robust framework for understanding cellular heterogeneity, ChrisCell-Graph holds transformative potential for the discovery of cell states and their specific functions that were previously challenging to characterize.

We calculated the Mean Square Distance [83] between the vectors corresponding to different states and visualized this information in the distance confusion matrix presented in Supplementary Fig.6a. This visualization supports our earlier conclusions. Specifically, states 3302, 3430, and 3174 exhibit the closest distances among all states, while states 1474 and 1490 also display relatively close distances to one another. Notably, the distance between states 3174 and 3302 is closer than that between either of these states and 3430, corroborating the findings illustrated in Fig. 7c. We posit that this characterization provided by ChrisCell-Graph offers a more detailed and precise quantization method for measuring the distances between cell states.

2.10 Gene Regulatory Network Analysis with ChrisCell-Graph
Studying gene-gene interactions is a critical step in understanding GRN [26] . GRNs are complex networks in which genes interact through regulatory elements, such as transcription factors, regulatory RNAs, and signaling pathways. These networks regulate gene expression levels and timing, thereby influencing cellular function and behavior. Investigating GRNs is essential for elucidating biological processes, understanding disease mechanisms, and developing therapeutic strategies.

We propose that ChrisCell-Graph serves not only as a valuable tool for studying gene interactions through its graph attention module but also visualizes these interactions at the cell state level via the ChrisCell SCD module. To validate this, we applied the ChrisCell-Graph model to the human pancreas cell dataset, focusing on Acinar cells.

Fig. 8a displays the gene-gene interactions predicted by the ChrisCell-Graph for an acinar cell. We encoded the genes from the acinar cell and obtained their corresponding codes by SCD module, as SCD module can also embed the gene embeddings into a code. Two genes are defined as connected if they share the same code. We selected highly expressed genes in acinar cells to visualize their interactions, as shown in Fig. 8a.

The gene regulatory network predicted by ChrisCell-Graph is biologically meaningful. For instance, ChrisCell-Graph predicts an interaction between the PRSS1 and CTRB1 genes. PRSS1 encodes trypsinogen, a major digestive enzyme in the pancreas, while CTRB1 encodes a trypsin inhibitor that regulates the activity of pancreatic digestive enzymes and prevents their premature activation in alveolar cells. The expression and function of these genes are interdependent; overexpression of PRSS1 may result in insufficient inhibition by CTRB1, thereby increasing the risk of cellular injury. The importance of these genes in the Gene Ontology database is rated at 0.5, indicating a significant biological relevance.

Another example is the interaction between CTRC and CFB genes. CTRC encodes a regulator of chymotrypsinogen, a pancreatic protease involved in the activation and regulation of digestive enzymes. CFB encodes complement factor B, which plays a role in the immune response and may influence the inflammatory response in the pancreas. In pathological states such as pancreatitis, CTRC expression may be modulated, while CFB activity can affect the extent of the inflammatory response. Therefore, the interaction between these genes may contribute to the response and repair mechanisms of the pancreas. Their importance in the Gene Ontology database is rated at 0.363, indicating a critical co-expression relationship. These examples demonstrate that the ChrisCell-Graph effectively captures the co-expression relationships of genes, serving as a valuable analytical tool for discovering gene regulatory networks in individual cells.

Fig. 8b-d illustrates the gene interactions across three primary states. Fig. 8b presents the overall gene-gene interactions within state 1475 cells. Notably, we found several genes that exhibit interactions with numerous other genes. The first example is the TPT1 gene [84], also known as Tumor Protein 1 or TCTP (Translationally Controlled Tumor Protein). TPT1 encodes a widely expressed intracellular protein involved in various biological processes, including cell proliferation, survival, and stress response. It interacts with genes that regulate the cell cycle and influence proliferation rates [84]. TPT1 is connected to ribosomal protein family genes, such as RPL31, RPLP0, and RPS8, which are essential components of the ribosome and play critical roles in protein synthesis, displaying high attention weights [84]. Furthermore, TPT1 exhibits a significant attention weight of 0.157 in relation to the H19 ENSG00000130600 [85] gene. This non-coding RNA gene, located on human chromosome 11, is crucial for embryonic development and cellular proliferation and is hypothesized to be part of an imprinted gene [85]. Thus, H19 ENSG00000130600 and TPT1 are interconnected in processes related to cell proliferation, tumorigenesis, and stress response.

The second example focuses on the DLK1 gene [86]. DLK1 (Delta-like 1) encodes a cell surface protein that is crucial for embryonic development, stem cell fate determination, and cell differentiation. It is closely associated with the Notch signaling pathway, growth factors, transcription factors, and tumor-related genes [86]. As shown in Fig. 8b, DLK1 exhibits a high attention weight of 0.133 in relation to the CLU gene. The CLU gene (Clusterin) encodes a widely expressed glycoprotein that primarily mediates cellular stress responses, apoptosis, and immune responses [87]. CLU plays a significant role in various physiological and pathological processes, including neuroprotection and tumorigenesis [87]. Both CLU and DLK1 are implicated in the cellular response to environmental stress and are aberrantly expressed in numerous cancers, linking them to tumorigenesis and progression. Their interaction may involve the regulation of cell proliferation and apoptosis, with DLK1 potentially modulating the expression or function of CLU via the Notch signaling pathway, suggesting a mutually dependent relationship in cellular processes. Consequently, ChrisCell-Graph assigns a high attention weight to the interaction between these two genes. Additionally, DLK1 is connected to several ribosomal protein family genes, including RPS23, RPL11, and RPLP0, which are integral to cell development and differentiation [86].

The final example pertains to the intrinsic interactions among ribosomal protein family genes. For instance, genes such as RPS23, RPL8, RPL13, and RPLP0 are linked within the ChrisCell-Graph framework and exhibit high attention weights. This suggests that these ribosomal protein genes work together to form ribosomal complexes and influence each other’s expression, leading to increased attention weights.

Fig. 8c illustrates the gene-gene interactions in state 3826 cells. We observe that the gene-gene interactions in state 3826 cells do not significantly differ from those in state 1475 cells, suggesting that the gene-gene interactions within specific cell types remain relatively stable. However, we note that ribosomal protein family genes, such as RPS23 and RPL13, demonstrate more pronounced interactions with other genes in this state. This phenomenon is likely linked to the active differentiation or developmental processes characterizing state 3826 cells.

Fig. 8d presents the gene-gene interactions in state 820 maintain a similar pattern to those observed in the other two states. This further supports our earlier conclusion that, despite differences in cell states, the gene-gene interactions within a specific cell type adhere to gene regulatory network principles [88].

The gene regulatory network analysis conducted using ChrisCell-Graph reveals its capability to capture gene-level knowledge and connections, identifying crucial regulatory genes such as DLK1 and TPT1. Moreover, ChrisCell-Graph provides interaction analysis results at the cell state level, allowing users to explore not only state-gene interactions but also to uncover gene-gene interactions. This dual capability further enhances our understanding of cellular mechanisms.

2.11 Multi-modal analysis with ChrisCell-Graph
For cell state identification, manual annotations such as cell types and tissue characteristics are critical informational factors that contribute to the construction of cell states. However, prior research has predominantly concentrated on analyzing scRNA-seq or scATAC-seq data, often neglecting the influence of these properties. To address this gap, we introduce a scalable approach that employs a language tokenizer which introduced by BioGPT [89]. It include 42,383 tokens and can convert the text of single cell properties description into one or a seires of tokens, effectively embedding these properties within a unified graph framework. This enables a comprehensive examination of how these annotated properties influence cell states.

To evaluate the significance of these properties, we validated ChrisCell-Graph using a sampled pretraining dataset comprising 10,000 cells, categorized into 167 cell types, sourced from 30 different tissues and involving 465 donors. This dataset enables us to consider various properties, including tissue type, disease state, age, and sex, and to investigate their respective influences on cell states.

Fig. 9a-d depict the distribution of attention weights for Tissue, Disease, Age, and Sex across nine distinct states. Our analysis reveals that tissue annotations are the most influential factors among these properties, with attention weights ranging from 0 to 0.25 and an average attention weight of 0.12 across the states. Additionally, Fig. 9e illustrates a strong correlation between cell distribution and tissue types, indicating that cells sharing the same tissue tend to cluster together in specific areas of the visualization.

Following tissue annotations, disease annotations represent the second most significant factor influencing cell states. While one might intuitively expect disease to play a major role, the predominance of normal cells within our pretraining dataset (as depicted in Fig. 9f) diminishes the perceived importance of disease, resulting in average attention weights that span from 0 to 0.14. Furthermore, the attention weights associated with age and sex are comparatively lower, averaging only 0.03 and 0.01, respectively. Fig. 9g and h provide UMAP visualizations of ChrisCell-Graph representations labeled by Age and Sex, revealing a relatively random distribution. This observation further substantiates the minimal contribution of age and sex to the construction of cell states.

In summary, ChrisCell-Graph offers a valuable tool for biological researchers, enabling the analysis of multimodal properties and facilitating the identification of crucial factors that contribute to cell states. By incorporating these annotations, ChrisCell-Graph enhances our understanding of the complexities underlying cellular behavior and state identification.

Methods
4.1 Preprocess
To preprocess gene expression and scATAC data, we begin by applying normalization and log1p transformation. These steps are essential for addressing scale disparities and mitigating the impact of missing values within the expression matrix. For the property data, which are typically represented as textual annotations, we utilize the BioGPT tokenizer [89] to convert these textual properties into a sequence of integers. This method enables us to leverage advanced language modeling techniques, allowing us to effectively represent cellular properties as embeddings.

4.2 ChrisCell Model Architecture
4.2.1 Input Embeddings
We utilize an embedding module, which consists of three components: gene expression, gene position (gene id), and gene ontology embeddings, to embed the normalized gene expression values. Given a series of gene expression values E = {e ,e ,e ...e } where 1 2 3 n
e i ∈ R, and a series of gene position indexes P input = {p 1,p 2,p 3...p n} where p i ∈ N.
We utilize a gene expression embedding module f and a gene position embedding module f to embed these two inputs:
E input =f E (e 1,e 2,e 3...e n) where E input ∈Rn×d (1)
P input =f P (p 1,p 2,p 3...p n) where P input ∈Rn×d (2)
where d is the embedding dimension of the model input.
Besides, we utilize the prior gene ontology knowledge to construct a symmetric importance matrix A, where A ∈Rn,n and A i,j represents the importance of gene i for gene j. We utilize a learnable matrix W A ∈Rn,d to embed the gene ontology matrix:
A input =AW A where A input ∈Rn,d (3)
Notably, we make the gene ontology matrix A learnable, meaning the gene ontology updates during training. Subsequently, we add these three embeddings as the model input:
X input =E input +P input +A input (4)

4.2.2 Encoder
To refine the gene features to a greater extent, we utilize multiple transformer layers to incorporate the information of expressed genes. Each transformer layer consists of two modules: the Self-Attention Mechanism and the Feed-Forward Neural Network.
Self-Attention Mechanism Given the encoder input X input = {x 1, x 2, x 3...x n} where X input ∈ Rn×d, the self-attention operation computes attention scores between all pairs of elements using query (Q), key (K), and value (V ) matrices. These matrices are linear transformations of the input embeddings:
Q=X input W Q, K =X input W K, V =X input W V . (5)
The attention weights are calculated as:
Xˆ =softmax( QKT / \sqrt{D} )V. (6)
This mechanism allows each gene to capture information from other genes.
Feed-Forward Neural Network After getting the output of Self-Attention Mechanism, we apply a Feed-Forward Neural Network F, which is a multi-layer perception function. What’s more, we add a residue connection:
X output =F(Xˆ)+X input . (7)
We use the X output as the next transformer layer’s input till get the final encoder output. Finally, we use a max polling to obtain the cell embeddings:
X cell =MaxPooling(X output ) (8)
The max pooling can capture the significant feature among all genes for each dimension.

4.2.3 Single-Cell Discretization Module
Following the Foldtoken [90], we adopt the Vector Quantization method called Soft-CVQ, which enables the transformation of continuous embeddings into discrete codes. Unlike traditional approaches that directly map embeddings to tokens, SoftCVQ employs a more nuanced mechanism to capture the relationships between embeddings. The process involves two key steps: (1) mapping pre-defined binary vectors (b j) to continuous token embeddings (v j) and (2) performing soft alignment between these token embeddings and latent embeddings (h s).
Given a decimal integer z and a codebook size m, the integer z is represented as a binary vector b j with a length of log 2(m). For example, if m=4, the binary vectors are b 1=[0,0], b 2=[0,1], b 3=[1,0], and b 4=[1,1]. The quantization operation converts continuous embeddings h s into discrete tokens z through the following process:
a sj = \sum j e = x 1 p e ( x f p s T ( v f j s T / v T j ) /T) (9)
z =argmax a sj
where v j represents the j-th token embedding, which also serves as the de-quantization component. During training, the binary vector b s is derived using bitwise operations:
b s =Bit(z,log 2(m)) (10)
fˆ s =ConditionNet(b s)
The temperature parameter T in Equation 9 controls the softness of the attention mechanism. When T is large, the attention weights approximate uniformity, whereas smaller T values result in weights closer to one-hot distributions. To stabilize training, T is gradually reduced from 1.0 to 10−8. The conditional network, ConditionNet, is a multi-layer perceptron (MLP) that maps high-dimensional binary vectors to latent embeddings of dimension d. For a codebook size of 216, this MLP transforms 65536 16-dimensional binary vectors into 65536 d-dimensional embeddings.
The core process of the aforementioned discretization encoding involves mapping continuous biomolecular data to a finite discrete binary encoding space, achieving efficient feature compression. This universal single-cell compilation method based on binary discrete encoding ensures that the model provides uniformly granular discrete representations for different genes and single-cell batches. It fully captures key biomolecular features, enabling standardized representation of heterogeneous data. Through a unified codebook, heterogeneous data can be compared, analyzed, and modeled within the same space, enabling seamless integration of data across different batches. We provide the discussion about how SCD module resolve the batch effect in Supplementary Section 11.

4.2.4 Decoder
After getting the encoder output, we concatenate the non-expressed, non-sampled gene embeddings, the encoder output, and the SCD embeddings as the decoder input.
X dec =[X ne ,X ns ,X SCD gene ,X SCD cell ]. (11)
Given that the decoder input consists of nearly 20,000 genes—too many for a standard Transformer to handle—we opted to use a Performer as the decoder. The Performer is a Transformer-based variant designed to enhance the computational efficiency of the self-attention mechanism. It achieves this goal by utilizing kernelized attention, which reduces computational complexity and enhances efficiency when handling long sequences. Given the decoder input X dec , a performer layer perform:
Q=X dec W Q, K =X dec W K, V =X dec W V . (12)
Xˆ dec = \sum n j=1 \phi(Q,K)V (13)
where \phi is a kernel function, which is used for approximating the attention matrix.
Then, we also perform a Feed-Forward Neural Network F
X decoutput =F(Xˆ dec ), (14)
After passing through six performer layers, we get the decoder outputs.

4.2.5 ZINB predictor
To predict the Zero-Inflated Negative Binomial, we deploy a ZINB predictor. Given input features X decoutput , the model output can be represented as:
Z ∼Bernoulli(π) (15)
where π is the probability of the zero-inflated part, typically predicted by logistic regression:
π =σ(W π X decoutput +b π ) (16)
Next, the parameters of the negative binomial distribution r and p are predicted by another regression model:
Y ∼NB(r,p) (17)
r =f r (X), p=f p (X decoutput ) (18)
where the NB denotes the negative binomial distribution. The final ZINB output can be expressed as:
Y = 0 with probability π
NB(r,p) with probability (1−π) (19)

4.3 ZINB loss funtion
We utilize the ZINB distribution predicted by the ZINB predictor and the ground truth to calculate the ZINB loss to train our models. Notably, we select the genes that are not sampled for the encoder input and a number of non-expressed genes, where the number of non-expressed genes is less than or equal to the non-sampled genes, to avoid the unbalanced distribution. Given observed data Y and predicted parameters π,r,p, the ZINB loss can be expressed as:
L(Y,π,r,p)=− \sum i=1 N ( Y i log(1− π i )+(1−Z i ) log(π i )+log \Gamma(Y i +r) / \Gamma(r)\Gamma(Y i +1) (20)
+ r log( p / p+1 ) + Y i log( 1 / p+1 ) ) (21)
Where: - Y i is the observed count for the i-th observation. - Z i is an indicator variable that is 1 if Y i = 0 (zero-inflated), and 0 otherwise. - Γ is the gamma function.
How ZINB predictor and loss function overcome the data sparsity? As prior works supported [16, 91], ZINB predictor effectively addresses data sparsity by combining a zero-inflation model with a negative binomial distribution. It first models additional zeros through the zero-inflation component, enabling the model to capture the high proportion of zeros in real data rather than treating them solely as outcomes of the negative binomial distribution. The negative binomial component models the variability and discreteness of non-zero count data, providing a more accurate description of sparse data. This combination makes ZINB particularly well-suited for handling highly sparse data in fields such as gene expression, thereby enhancing the model’s fitting capability and predictive accuracy.

4.4 ChrisCell-Graph model architecture
4.4.1 Graph construction
To identify the neighboring cells, we compute the mean squared distance between the gene expression profiles of a target cell and other cells within the same batch. We then select the top 15 closest cells based on this distance as its neighbors. To optimize computational efficiency, we first reduce the dimensionality of our analysis by focusing on the top 500 highly expressed genes, determined through covariance analysis of the expression data.
Furthermore, we conceptualize each gene and property as nodes within a network. For each cell, we isolate the top 200 genes based on expression levels, alongside chromatin peaks characterized by the highest 1000 peak values. This strategy establishes connections between the cellular representations and the most expressive genes and chromatin regions, thereby enhancing our model’s ability to capture relevant biological interactions.

4.4.2 Graph attention module
We introduce the graph attentional layer [54] designed for modeling relationships in graph-structured data, which is a fundamental component in various graph-based machine learning models. This layer is versatile and can be integrated into different Graph Attention Network (GAT) architectures, offering flexibility without assuming any specific attention mechanism.
Input and Transformation: The layer processes a set of node features h = {h 1,h 2,...,h N}, where each node i has features h i ∈RF. The transformation begins with a shared linear layer using a weight matrix W ∈RF′×F, projecting each node’s features to a new space Wh i ∈RF′.
Graph Attention Mechanism: The core of the layer is the attention mechanism, which computes coefficients e ij representing the importance of node j’s features for node i:
e ij =a(Wh i,Wh j) (22)
These coefficients are normalized using the softmax function to ensure they sum to 1 for each node (i):
α ij =softmax j (e ij )= exp(e ij ) / \sum k∈Ni exp(e ik ) (23)
where N i denotes node i’s neighborhood.
Feature Transformation: The new features for each node are computed as a weighted sum of its neighbors’ transformed features:
h i =σ( \sum j∈Ni α ij Wh j) (24)
Here, σ introduces nonlinearity, enhancing the model’s ability to capture complex patterns.
Multi-Head Attention: To enhance stability and model capacity, the layer employs multi-head attention. (K) independent attention mechanisms compute their own coefficients αk ij and apply their own transformations with weight matrices Wk. The outputs are concatenated:
h i = ∥ k=1 K σ( \sum j∈Ni α i k j Wkh j) (25)
This increases the model’s capacity by capturing diverse aspects of the data.
Final Layer Considerations: When applied to the final layer of a network, concatenation may not be optimal. Instead, averaging across heads is preferred:
h′ i =σ( 1 / K \sum k=1 K \sum j∈Ni α i k j Wkh j) (26)
This approach prevents over-reliance on any single head, improving generalization.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This manuscript presents ChrisCell, a single-cell foundation model with an integrated discretization module (SCD) and a graph-based extension (ChrisCell-Graph). While the scale of pretraining (68M cells) and claimed performance improvements are notable, the work suffers from fundamental issues in conceptual clarity, methodological justification, and biological interpretability. The presentation obscures technical contributions, and key components lack biological or computational motivation. These flaws significantly undermine the manuscript's scientific rigor and reproducibility.\n1. The proliferation of undefined terminology creates substantial barriers to understanding. Novel concepts like \"Single-Cell Discretization (SCD)\" are introduced without clear justification. The manuscript fails to establish why continuous embeddings require discretization, what specific limitations SCD resolves that existing methods cannot address, or how quantization biologically correlates with cellular states. Without ablation studies demonstrating performance degradation when removing SCD, its necessity remains unsubstantiated. This pattern extends to terms like \"cell code\" and \"data-driven identifier,\" which are neither mathematically defined nor biologically contextualized.\n2. The arbitrary numerical labeling of cell states (e.g., state 3430/3302) disregards established biological conventions. Cell states are intrinsically linked to functional, morphological, or molecular phenotypes—not abstract integers. The proposed \"state genes\" further conflates nomenclature; genes are experimentally established as markers of states via orthogonal validation (e.g., proteomics, spatial imaging). The absence of integration with modalities like Cell Painting or perturbational validation reduces these findings to computational artifacts. If states represent biological reality, they must align with measurable phenotypic axes.\n3. The introduction's three limitations (sparsity, heterogeneity, interpretability) are inadequately linked to the model's architecture. The manuscript does not specify how SCD mitigates heterogeneity versus standard batch correction, nor does it clarify which modules target sparsity versus interpretability. Claims about interpretability are circular: SCD's \"cell codes\" are deemed interpretable because they are discrete, yet no experiment demonstrates biological insight derived solely from these codes. The disconnect between stated limitations and architectural solutions suggests post-hoc justification.\n4. Methodological descriptions are critically deficient. The transformation of gene embeddings to cell embeddings via \"pooling\" lacks mathematical specification (e.g., max/mean/sum pooling?) and biological rationale—genes are not positionally invariant like pixels, so pooling operations require justification. Gene position indexes (Pinput) are undefined: do these represent genomic coordinates, canonical gene ordering, or random indices? The \"cell code\" generation (Eq 4-5) operates as a black box, with no discussion on how binary vectors map to biological meaning. Such omissions prevent replication or meaningful evaluation.\n5. The model's design lacks coherent motivation. Figure 1a depicts components (e.g., gene ontology embedding, SCD module) without explaining their necessity. Why is a Performer decoder used instead of a standard Transformer? How does discretizing cell embeddings aid downstream graph construction? The workflow appears as a patchwork of techniques rather than a principled solution to defined biological challenges. This disconnection between motivation and implementation extends to ChrisCell-Graph, where graph attention weights are interpreted as biological importance without statistical or experimental validation.\n6. The claim of interpretability for ChrisCell-Graph is unsubstantiated. Attention weights in GNNs are not intrinsically interpretable as biological mechanisms—they require rigorous benchmarking against known interactions. The plugin-based nature of ChrisCell-Graph undermines the core model's novelty: if any foundation model's embeddings can generate similar graphs, ChrisCell's value diminishes. Crucially, no comparisons against scGPT/scFoundation embeddings processed by ChrisCell-Graph are provided, leaving the central \"interpretability\" claim unsupported.", "reviewer_2": "Wang, Tan et al. present ChrisCell, a novel single-cell foundation model. The authors report that ChrisCell achieves state-of-the-art performance across multiple downstream tasks, including clustering, cell type annotation, and prediction of drug and perturbation responses. Given the rapid expansion of single-cell foundation models as a key area of AI in systems biology, the manuscript addresses a timely and broadly relevant topic. Nonetheless, I have several questions regarding both the methodological details and the scope of the applications presented.\nMajor\n1. Ablation studies: Critical ablations are missing to find out what makes ChrisCell perform well. The model does differ a lot from the \"mainstream\" foundation models e.g. scGPT, NicheFormer as these are usually decoder-only models which represent the masked genes in the same context as the non-masked. To understand what are the critical components of ChrisCell, ablations are needed as described in more detail below:\n* SCD module\n* Less parameters: ChrisCell is huge (500M params) and at least transformer layers need to be ablated\n* ZINB decoder\n2. Clustering and cell annotation: while the authors show superior performance on the investigated datasets, some more detailed benchmarking (including other baseline models and single cell datasets) would be beneficial. Recent benchmarking efforts, like the ones from CZI (https://virtualcellmodels.cziscience.com/benchmark-task/cell-type-classification) can provide additional datasets and models for this.\n3. Drug sensitivity prediction: Importantly, the description and the citation to the used drug sensitivity datasets is missing. It is also not clear, which splits the authors used in the drug response classification problem. As SCAD performs very badly on this task, some additional baseline model would be needed. For the drug response prediction (IC50 regression) task the authors used random split of the train and test data. Random split, especially with IC50 data (Ovchinnikova et al. npj Precision Oncology 2024) can lead to irrealistic high performance metrics, and can mask real differences between models. To show the superior performance of ChrisCell, some more realistic train / test splits (new cell lines, new drugs etc.) would be required.\n4. Perturb-seq prediction: Recently several papers showed that simple baselines (like linear or random forest models, some examples in Ahlmann-Eltze et al. Nature Methods 2025), and even the simplest train mean model outperform foundation models for post-treatment RNAseq prediction task. The authors should also compare their results to these simple models.\n5. GRN analysis: The identified interactions (like PRSS1 and CTRB1) are probably not gene regulatory, but simple gene co-expression interactions, as none of the identified genes are transcription factors. The authors should compare their gene interaction network with some gene regulatory database, like CollecTRI (Muller-Dott et al. 2023 Nucleic Acids Research).\nMinor\n1. What is meant by \"setting the SCD cell embedding at the ”Count” token for the cell prompt\". The cell embedding seems to be \"just\" an embedding concatenated to the decoder context and not used in a distinguished way.\n2. Some discussion would be desired on how the \"mainstream\" cell foundation models are different from ChrisCell e.g. decoder-only, masked gene expression / gene pretraining task. It would be fruitful to understand the authors’ intuitions.\n3. ChrisCell embeddings were used as a replacement for GO graph of co-expression graph of GEARS?\n4. Based on my knowledge, the Adamson perturb-seq dataset contains only 100 different perturbations. In Figure 4E, on the left panel there are many more data points. What does each data point represent? A pseudobulked perturbation (as in the original GEARS evaluation), or individual single cells?\n5. The authors write: “Fig. 5b presents the overall performance of ChrisCell and other baseline methods across three perturb-seq datasets, evaluated using four metrics.” - in the drug (not perturb-seq) perturbation prediction section.", "reviewer_3": "Wang et al. present ChrisCell as a novel single-cell foundation model by combining VQ-VAE. The model was shown to outperform other large language models on several single-cell data analysis tasks. However, the comparative analysis is not comprehensive enough, as only a limited number of baseline methods are included. The introduction of ChrisCell-Graph as an interpreter is interesting, but its actual contribution to the overall interpretability of the model remains unclear. Many technical details of ChrisCell-Graph are not sufficiently explained and current results can’t be repeated or well evaluated.\nI have the following major concerns:\n1. Lack of Comprehensive Comparison with Existing Methods. The comparative study with scGPT and scFoundation is narrow. The author need include recent large language models, such as UCE[1], GeneCompass[2], and CellFM[3] for comparison.\n2. Insufficient Justification of Innovation. The core innovation is stated to be the integration of SCD. However, the improvements brought by SCD are unclear, and it is also not well demonstrated how SCD addresses the challenge of heterogeneous datasets. The authors should provide a comparison between the encoder-decoder and encoder-SCD-decoder architectures under similar parameter settings.\n3. Unclear Description and Application of ChrisCell-Graph. The manuscript provides a confusing analysis of ChrisCell and ChrisCell-Graph, making it difficult to understand the relationship and distinctions between the two modules. It is unclear which tasks use ChrisCell and which use ChrisCell-Graph. For example, Section 2.3.1 introduces ChrisCell-Graph, but the identification of marker genes only uses the embeddings from ChrisCell; Section 2.3.3 is titled \"Gene Regulatory Network Analysis with ChrisCell-Graph,\" but the content does not mention ChrisCell-Graph at all or explain how it is applied in this context.\n4. Please specify the tasks involving fine-tuning and those based solely on embedding extraction. Are the same settings used for the comparison with other foundation models? The comparisons may not be fair if your model used the finetuning while other models used only the embedding.\n5. The author may consider including more baseline methods of cell type prediction (See references https://genomebiology.biomedcentral.com/articles/10.1186/s13059-021-02480-2 and https://genomebiology.biomedcentral.com/articles/10.1186/s13059-019-1795-z). In addition, scBERT has shown excellent performance in cell type annotation. Does ChrisCell also achieve superior or at least comparable performance to scBERT in this task?\n6. 4.3. In the drug response classification task, SCAD is trained on labeled bulk RNA-seq data and unlabeled single cell data. It is not clear how ChrisCell is trained and evaluated in this setting.\n7. In the drug perturbation prediction task, please clarify how the data are split for training and testing. Additionally, it would be beneficial to include comparisons with several non-foundational models, such as scVIDR[4] or CellOT[5].\n8. The meaning of Supplementary Fig. 3a is unclear. Please clarify what this figure represents and explain how it supports the claim of validating the effectiveness of ChrisCell-Graph’s attention weights.\n9. The process for generating Figure 7d is unclear and requires further explanation.\n10. The term \"language tokenizer\" is introduced but not sufficiently explained in the manuscript.\n11. In Figure 9a-d, the meaning of the x-axis is unclear. Additionally, please clarify how the attention score is calculated and what entities are involved in this calculation.\n12. In section 2.3.5, the procedure for fine-tuning on ATAC-seq data should be described in detail.\nOther issues:\n13. Most existing large models use MSE as the reconstruction loss, while the authors use ZINB loss. What is the specific benefit of ZINB over MSE? It would be better to provide an ablation experiment comparing the two losses to clarify the impact.\n14. There is an unexpected red wavy line appearing in Supplementary Figure 1.\n15. In Section 2.3.2 (last paragraph) and Section 2.3.3 (third paragraph), figure references are incorrect.\n16. Some biological insights require literature support, such as the marker genes shown in Figure 6c.\n[1] Rosen Y, Roohani Y, Agarwal A, et al. Universal cell embeddings: A foundation model for cell biology[J]. bioRxiv, 2023: 2023.11. 28.568918.\n[2] Yang X, Liu G, Feng G, et al. GeneCompass: deciphering universal gene regulatory mechanisms with a knowledge-informed cross-species foundation model[J]. Cell Research, 2024, 34(12): 830-845.\n[3] Zeng Y, Xie J, Shangguan N, et al. CellFM: a large-scale foundation model pre-trained on transcriptomics of 100 million human cells[J]. Nature Communications, 2025, 16(1): 4679.\n[4] Kana O, Nault R, Filipovic D, et al. Generative modeling of single-cell gene expression for dose-dependent chemical perturbations[J]. Patterns, 2023, 4(8).\n[5] Bunne C, Stark S G, Gut G, et al. Learning single-cell perturbation responses using neural optimal transport[J]. Nature methods, 2023, 20(11): 1759-1768.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "NCOMMS-25-59714-T response\nWe thank all reviewers for the insightful suggestions and comments to help us improve\nthe manuscript. Below are our point-by-point responses to address all the concerns\nraised by reviewers. All the changes made in the main manuscript are highlighted in\nyellow.\nOverview of the Revisions:\n1. We provided a comprehensive ablation study to investigate the contribution of\neach module\nIn response to Questions 1, 5, and 6 from Reviewer 1, Question 1 from Reviewer 2,\nand Question 2 and 13 from Reviewer 3, we conducted a comprehensive ablation\nstudy to evaluate the individual contributions of each module in our proposed\nmodels. The study includes an analysis of the SCD module, GO embedding module,\nZINB decoder, ChrisCell models with 50M and 100M parameters, as well as the\nChrisCell-Graph model that integrates embeddings from multiple foundation\nmodels.\n2. We conducted the in silico perturbation experiments to validate the cell states\nand state genes provided by ChrisCell and ChrisCell-Graph.\nIn response to Question 2 and 5 from Reviewer 1, we employed CellOracle to\nconduct in silico perturbation predictions. The results demonstrate that cells within\nthe same cluster exhibit significant biological commonality, and that state-specific\ngenes substantially influence cell state construction and transition.\n3. We applied ChrisCell on more cell clustering and cell type prediction\nbenchmarks to further confirm the generalizability of ChrisCell.\nIn response to Question 2 from Reviewer 2 and Question 5 from Reviewer 3, we\nevaluated ChrisCell on the cz-benchmark, which comprises five cell clustering and\nfive cell type prediction tasks. Notably, ChrisCell consistently excelled across all\nthese tasks.\n4. We added more baselines on each downstream task to further demonstrate the\nadvance of ChrisCell.\n`\nIn response to Questions 2, 3, and 4 from Reviewer 2 and Questions 1 and 5 from\nReviewer 3, we have expanded the set of baseline methods across each downstream\ntask.\nFor cell clustering, cell type annotation, and property prediction tasks, we\nadditionally included GeneCompass, CellFM, UCE, and scBERT. We also\nincorporated scVI, TranscriptFormer, Geneformer, and AI.DO on the cz-benchmark.\nFor gene perturbation tasks, we further added Linear, GeneCompass, and CellFM\nacross three datasets. For drug response prediction tasks, we supplemented the\ncomparisons with Linear and Random Forest algorithms.For drug perturbation tasks,\nwe additionally evaluated scVIDR and CellOT. Furthermore, we demonstrated that\nChrisCell continuously achieve superior performance.\n5. We compared ChrisCell-Graph with a gene regulatory database of CollecTRI.\nIn response to Questions 5 from Reviwer 2, we compared ChrisCell-Graph with the\ngene regulatory database of CollecTRI and show that ChrisCell-Graph can provided\nsimilar but also insightful gene regulatory network, demonstrating that ChrisCell-\nGraph is also capable of capturing the TF regulatory pattern.\n6. We added detailed describtion and discussion of our methods to improve the\ncomprehensibility and reproducibility.\nIn response to Question 3-5 from Reviewer 1 and Question 2-4 and 8-12 from\nReviwer 3, we further included more extensive description and discussion about our\nmethods, including how SCD module overcome the data heterogeneity, ZINB solve\nthe data sparsity, and the contribution of each module, and some details about the\nexperimental settings and method describtion.\n7. We made a comprehensive comparison of the ChrisCell architecture with other\nfoundation models.\nIn response to minor concern 2 from Reviewer 2, we conduct a comprehensive\narchitecture comparison of ChrisCell with other foundation models, including both\ngenerative (scGPT) and masked (scBERT, scFoundation) Single-Cell Foundation\nModels, to show the uniqueness of ChrisCell.\n`\nRespond List to Comments:\nReviewer #1 (Remarks to the Author):\nThis manuscript presents ChrisCell, a single-cell foundation model with an integrated\ndiscretization module (SCD) and a graph-based extension (ChrisCell-Graph). While\nthe scale of pretraining (68M cells) and claimed performance improvements are\nnotable, the work suffers from fundamental issues in conceptual clarity, methodological justification, and biological interpretability. The presentation obscures technical contributions, and key components lack biological or computational motivation. These flaws significantly undermine the manuscript's scientific rigor and reproducibility.\n1. The proliferation of undefined terminology creates substantial barriers to\nunderstanding. Novel concepts like \"Single-Cell Discretization (SCD)\" are introduced\nwithout clear justification. The manuscript fails to establish why continuous\nembeddings require discretization, what specific limitations SCD resolves that\nexisting methods cannot address, or how quantization biologically correlates with\ncellular states. Without ablation studies demonstrating performance degradation\nwhen removing SCD, its necessity remains unsubstantiated. This pattern extends to\nterms like \"cell code\" and \"data-driven identifier,\" which are neither mathematically\ndefined nor biologically contextualized.\nResponse: Thank you for highlighting this important issue!\nThe primary motivation for introducing the SCD module is to overcome the critical\nchallenge of data heterogeneity in large-scale single-cell datasets. Data heterogeneity\nis a well-recognized challenge across numerous database fields. In single-cell genomics,\nit arises predominantly from batch effects—technical variations introduced when data\nare generated across different sources, platforms, or experimental protocols. These\ndiscrepancies in processing methods and measurement scales result in systematic biases\n`\nthat obscure biological signals.\nThe reason of discretization is that when single-cell foundation models are pre-\ntrained on large-scale datasets comprising multiple batches, they may learn to exploit\nbatch-specific artifacts rather than biologically meaningful patterns. This leads to\ninformation leakage: models can achieve high performance by recognizing batch\nidentifiers, especially since cells within a batch often share cell types or expression\nprofiles. However, such models generalize poorly to unseen data containing new\nbatches, as they fail to adapt to technical variations absent during training. We\ncompared the batch correction capability of ChrisCell and an SCD-ablated variant using\nboth a training subset (10,000 cells from 452 donors) and a held-out test set. As shown\nin Fig.5 c–d, the model without SCD achieved an AvgBatch score of 18.33 on the\ntraining data, but only 9.86 on the test set—indicating significant overfitting to batch\neffects. This conclusion is supported by the ablation results in Supplementary Fig.1,\nwhere the same model exhibited rapid training convergence yet delivered the poorest\ndownstream performance across all tasks.\nThe SCD module addresses this problem by projecting continuous cell and gene\n`\nembeddings into a discrete, unified vector space through vector quantization.\nDiscretization forces the model to compress the embeddings, filtering out noisy, high-\nfrequency variations that are often associated with batch effects, while preserving the\ncore, low-frequency signals corresponding to fundamental biological states. Moreover,\nit maps cells from different batches, with different expression scales and technical\nartifacts, into a shared, universal codebook. This leads us to our terminology: \"Cell\nCode\": This is the discrete vector from the unified codebook that represents a cell after\nquantization. \"Data-Driven Identifier\": We use this term to describe the \"cell code\"\nbecause it is a standardized identifier learned from the data itself, representing a cell's\nbiological state in a way that is robust to the technical artifacts of its original batch.\nParallel insights can be drawn from protein language modeling. For instance,\nAlphaFold2’s predictions can be influenced by template-induced bias, leading to\nstructural inaccuracies. Methods such as Saprot[1] and ESM3[2] introduced vector\nquantization (VQ) to project amino acid embeddings into a discrete, unified space,\nthereby reducing superficial biases while preserving essential biological information.\nAs a result, the SCD embeddings capture fundamental cellular features without batch\neffects. Empirical evidence supports this mechanism: ChrisCell demonstrates\nconsistent AvgBatch scores across both training and test sets (Fig.5c-d), and maintains\nrobust performance in ablation studies (Fig.1), confirming the efficacy of SCD in\novercoming data heterogeneity.\nWe conducted the comprehensive ablation studies, with full results presented in\nSupplementary Section 1 and a detailed discussion in Supplementary Section 9.\n`\nNotably, we find that the model without the SCD module obviously poses the worst\nperformance on all datasets and metrics. However, Fig. 1d illustrates that the model\nwithout the SCD module poses the best training convergence among all models. This\nphenomenon provides strong evidence that the model without the SCD module suffers\nfrom data heterogeneity, and the SCD module greatly overcomes this issue.\n2. The arbitrary numerical labeling of cell states (e.g., state 3430/3302) disregards\nestablished biological conventions. Cell states are intrinsically linked to functional,\nmorphological, or molecular phenotypes—not abstract integers. The proposed \"state\ngenes\" further conflates nomenclature; genes are experimentally established as\nmarkers of states via orthogonal validation (e.g., proteomics, spatial imaging). The\nabsence of integration with modalities like Cell Painting or perturbational validation\nreduces these findings to computational artifacts. If states represent biological reality,\n`\nthey must align with measurable phenotypic axes.\nResponse: Thank you for highlighting this important issue! We supplemented the\nrelevant experimental validation of in silico gene perturbation in the Supplementary\nSection 6, as follows:\nTo validate the cell states and state-specific genes identified by ChrisCell-Graph,\n`\nwe performed in silico perturbation experiments using CellOracle [3], a widely adopted\ncomputational model for perturbation prediction. We selected single cells assigned to\nstate 3174 from the human brain dataset described in Section 2.3.2 of the main text,\nresulting in 1103 oligodendrocytes, and systematically knocked out a set of genes to\nexamine changes in transcriptional activity.\nFig.7a shows the average gene expression differences between unperturbed and\nperturbed cells from state 3174, following the knockout of selected genes. The top eight\ngenes with the largest expression differences are highlighted. Among these, knockout\nof ANK3 and PCDH9—both identified as state genes by ChrisCell-Graph—resulted in\nthe most substantial expression changes, suggesting that these genes play pivotal roles\nin state 3174.\nFig.7b displays the top 20 differentially expressed genes between control and\nperturbed cells within state 3174. Again, ANK3 and PCDH9 exhibited the most\npronounced knockout effects, further supporting their functional significance.\nMoreover, as illustrated in Fig.7c, knockout of ANK3 and PCDH9 induced state\ntransitions in 76% and 68% of the cells, indicating that cellular identity in\nstate 3174 is highly sensitive to perturbations in these two genes.\nAdditionally, as shown in Fig.7d, ANK3 expression in state 3174 cells was consistently\nhigh (ranging from 0.6 to 2.2) and stable, reinforcing the importance of ANK3 in\nmaintaining this cell state.\nIn Fig.7e-h, we also plot the in silico perturbation results for state 1474 cells, which\nshowcase similar results. The DOCK4 and PAG1 genes, provided by ChrisCell-Graph,\nare highlighted as crucial genes for cells in state 1474, further indicating that ChrisCell-\nGraph indeed captures biological insights at the cell state level and identifies key genes\nsignificantly impacting cell states.\n3. The introduction's three limitations (sparsity, heterogeneity, interpretability) are\ninadequately linked to the model's architecture. The manuscript does not specify how\nSCD mitigates heterogeneity versus standard batch correction, nor does it clarify\nwhich modules target sparsity versus interpretability. Claims about interpretability\nare circular: SCD's \"cell codes\" are deemed interpretable because they are discrete,\n`\nyet no experiment demonstrates biological insight derived solely from these codes.\nThe disconnect between stated limitations and architectural solutions suggests post-hoc\njustification.\nResponse: Thank you for highlighting this important issue!\nRegarding to the sparsity, we mentioned that we utilize the ZINB predictor and\nloss function to solve the data sparsity in the Introduction Section:“At the gene level,\nwe implement a cell-gene-prompt (CGP)-based Zero-Inflated Negative Binomial (ZINB)\nprediction task, which significantly enhances cell representation and captures gene co-\nexpression patterns, and addresses the high sparsity of the data.”. And we have added\nadditional content to explain why ZINB predictor can resole the data sparsity in Section\n4.4 of the revised manuscript, as follows:\nAs prior works supported [4,5], ZINB predictor effectively addresses data sparsity\nby combining a zero-inflation model with a negative binomial distribution. It first\nmodels additional zeros through the zero-inflation component, enabling the model to\ncapture the high proportion of zeros in real data rather than treating them solely as\noutcomes of the negative binomial distribution. The negative binomial component\nmodels the variability and discreteness of non-zero count data, providing a more\naccurate description of sparse data. This combination makes ZINB particularly well-\nsuited for handling highly sparse data in fields such as gene expression, thereby\nenhancing the model's fitting capability and predictive accuracy.\nRegarding to the data heterogeneity, we extend the explanation of how SCD\nmodule overcome the data heterogeneity in Supplementary Section 9, and we also have\ncomprehensively discussed it refer to the question 1.\nRegarding to the data in interpretability of the cell codes, in our manuscript, we\nposit that the cell codes generated by ChrisCell enable the identification of cellular\nfunctions and activities at the single-cell level, going beyond coarse cell-type\ncategorization. To support this claim, we performed substantial experiments in Section\n2.3, integrating ChrisCell’s embeddings and cell codes with ChrisCell-Graph to analyze\nand interpret biologically meaningful insights from real-world data. These results\ndemonstrate that the cell codes produced by ChrisCell are not only computationally\n`\ndiscriminative but also reflect relevant biological characteristics.\n4. Methodological descriptions are critically deficient. The transformation of gene\nembeddings to cell embeddings via \"pooling\" lacks mathematical specification (e.g.,\nmax/mean/sum pooling?) and biological rationale—genes are not positionally\ninvariant like pixels, so pooling operations require justification. Gene position indexes\n(Pinput) are undefined: do these represent genomic coordinates, canonical gene\nordering, or random indices? The \"cell code\" generation (Eq 4-5) operates as a black\nbox, with no discussion on how binary vectors map to biological meaning. Such\nomissions prevent replication or meaningful evaluation.\nResponse: Thank you for highlighting this important deficiency.\nRegarding the pooling method, we employed max pooling to generate cell\nembeddings. This approach effectively captures the most salient features across all\ngenes for each dimension of the representation.\nConcerning the gene position indexes, these correspond to specific gene\nidentifiers (such as CLU and DLK1). To prevent confusion, we have added clarification\nin the revised manuscript.\nAs for the SCD module, we have expanded the discussion in Section 4.2.3. The\ncore of the discretization encoding process involves mapping continuous biomolecular\ndata into a finite binary encoding space, enabling efficient feature compression. This\nuniversal single-cell compilation method, grounded in binary discrete encoding,\nensures that the model produces consistently granular discrete representations across\ndiverse genes and single-cell batches. It effectively captures essential biomolecular\ncharacteristics and facilitates standardized representation of heterogeneous data. By\nleveraging a unified codebook, heterogeneous datasets can be compared, analyzed, and\nmodeled within a shared space, thereby supporting seamless integration across different\nbatches.\n5. The model's design lacks coherent motivation. Figure 1a depicts components (e.g.,\ngene ontology embedding, SCD module) without explaining their necessity. Why is a\nPerformer decoder used instead of a standard Transformer? How does discretizing\ncell embeddings aid downstream graph construction? The workflow appears as a\n`\npatchwork of techniques rather than a principled solution to defined biological\nchallenges. This disconnection between motivation and implementation extends to\nChrisCell-Graph, where graph attention weights are interpreted as biological\nimportance without statistical or experimental validation.\nResponse: Thank you for pointing this issue! For gene ontology embedding module,\nwe mentioned in the Section 2.1:“The Gene Ontology embedding module incorporates\nprior knowledge from Gene Ontology to model gene co-expression”. Furthermore, we\nconduct comprehensive ablation studies to investigate the contribution of gene ontology\nembedding module in Supplementary Section 1, as follows:\nIn Fig.1, the model without the Gene Ontology (GO) embedding module exhibits\nthe smallest performance degradation compared to the baseline, indicating that the GO\n`\nembedding module is the least impactful module among other modules. However, we\nobserve that the performance of the model without the GO embedding module degrades\nmore significant in the gene perturbation tasks. For example, the PCC Delta of the\nmodel without GO embedding on three datasets is 0.601, 0.432, and 0.412, which is\ncomparable to or less than the model without the ZINB loss function(0.610, 0.431,\n0.409), while in other tasks, the model without the GO embedding module is\nsignificantly better than the model without ZINB loss. We attribute this phenomenon to\nthe prior gene ontology knowledge that helps model better capture the gene-gene\ninteraction, leading to better performance of gene perturbation prediction.\nFor SCD module, we have clarified the contribution and the necessity of SCD\nmodule in Supplementary Section 1 and Section 9 in the revised manuscript, and can\nrefer to the question 1.\nFor the Performer decoder, we mentioned in Section 4.2.4 of our manuscript:“It\nachieves this goal by utilizing kernelized attention, which reduces computational\ncomplexity and enhances efficiency when handling long sequences.” Given that the\ndecoder input consists of nearly 20,000 genes—too many for a standard Transformer\nto handle—we opted to use a Performer as the decoder. This clarification has been\nincluded in the revised manuscript in Section 4.2.4.\"\nFor the experimental validation of ChrisCell-Graph, we conducted the in silico\nperturbation experiments using CellOracle to verify the validity of attention weights\nprovided by ChrisCell-Graph in the Supplementary Section 6 in the revised manuscript\nand can refer to the question 2.\n6. The claim of interpretability for ChrisCell-Graph is unsubstantiated. Attention\nweights in GNNs are not intrinsically interpretable as biological mechanisms—they\nrequire rigorous benchmarking against known interactions. The plugin-based nature\nof ChrisCell-Graph undermines the core model's novelty: if any foundation model's\nembeddings can generate similar graphs, ChrisCell's value diminishes. Crucially, no\ncomparisons against scGPT/scFoundation embeddings processed by ChrisCell-Graph\nare provided, leaving the central \"interpretability\" claim unsupported.\nResponse: Thank you for raising this problem! We conduct the comparison that\n`\ncombining three foundation models with ChrisCell in the Supplementary Section 7,\nas follows:\nTo assess the capability of ChrisCell-Graph in leveraging embeddings produced by\nChrisCell, we conducted comparative experiments using embeddings from scGPT and\nscFoundation to reconstruct knowledge graphs. These graphs were generated and\nevaluated under the same framework applied to human brain single-cell data in Section\n2.3.2, with resulting sub-graphs illustrated in Fig.7.\nFig.7a presents the knowledge graph derived from ChrisCell embeddings using\nChrisCell-Graph. The model successfully associates cell types with functionally\nrelevant genes through high-edge-weight connections—for instance, linking\noligodendrocytes to ANK1 and PCDH9, and neurons to ANKS1B and RBFOX1.\nFurthermore, ChrisCell-Graph not only clusters cells of the same type and captures\nproximal relationships between neurons and oligodendrocytes but also clearly separates\ndistinct cell sub-graphs, such as those of oligodendrocytes and central nervous system\n`\nmacrophages. Additionally, the model resolves finer-grained cellular states: for\nexample, oligodendrocytes in state 3174 show stronger association with ANK3 and\nPCDH9, while those in state 3302 exhibit higher edge weights for FTX and ANK3, and\ngreater proximity to neuronal clusters—consistent with findings presented in Section\n2.3.2.\nAs shown in Fig.7b and c, knowledge graphs constructed from scGPT and\nscFoundation embeddings using ChrisCell-Graph also connect cells of the same type to\nrelevant genes, confirming the plug-and-play adaptability of our graph inference\nframework. However, edge weights in these graphs are consistently low (0.05–0.11),\nfailing to emphasize biologically meaningful relationships. Moreover, both graphs\nlack clear separation between cell-type-specific sub-graphs, resulting in\nambiguous structural representation. In contrast, ChrisCell produces a well-\nseparated, interpretable graph with salient biological connections. These results indicate\nthat ChrisCell generates more discriminative and informative representations than\nexisting foundation models, and affords enhanced resolution for cell-state-level\nanalysis, underscoring its superior interpretability.\nReference\n[1] Su, Jin, et al. \"Saprot: Protein language modeling with structure-aware vocabulary.\" BioRxiv\n(2023): 2023-10.\n[2] Hayes, Thomas, et al. \"Simulating 500 million years of evolution with a language model.\"\nScience 387.6736 (2025): 850-858.\n[3] Kamimoto, Kenji, Blerta Stringa, Christy M. Hoffmann, Kunal Jindal, Lilianna Solnica-Krezel,\nand Samantha A. Morris. \"Dissecting cell identity via network inference and in silico gene\nperturbation.\" Nature 614, no. 7949 (2023): 742-751.\n[4] Kalfon, J., Samaran, J., Peyre, G., Cantini, L.: scprint: pre-training on 50million cells allows\nrobust gene network predictions. Nature Communications 16(1), 3607 (2025)\n[5] Lopez, R., Regier, J., Cole, M.B., Jordan, M.I., Yosef, N.: Deep generative modeling for single-\ncell transcriptomics. Nature methods 15(12), 1053-1058 (2018)\n`\nReviewer #2 (Remarks to the Author):\nWang, Tan et al. present ChrisCell, a novel single-cell foundation model. The authors\nreport that ChrisCell achieves state-of-the-art performance across multiple\ndownstream tasks, including clustering, cell type annotation, and prediction of drug\nand perturbation responses. Given the rapid expansion of single-cell foundation models as a key area of AI in systems biology, the manuscript addresses a timely and\nbroadly relevant topic. Nonetheless, I have several questions regarding both the\nmethodological details and the scope of the applications presented.\nMajor\n1. Ablation studies: Critical ablations are missing to find out what makes ChrisCell\nperform well. The model does differ a lot from the \"mainstream\" foundation models\ne.g. scGPT, NicheFormer as these are usually decoder-only models which represent\nthe masked genes in the same context as the non-masked. To understand what are the\ncritical components of ChrisCell, ablations are needed as described in more detail\nbelow:\n* SCD module\n* Less parameters: ChrisCell is huge (500M params) and at least transformer layers\nneed to be ablated\n* ZINB decoder\nResponse: Thank you for highlighting this important issue! We supplemented this\ncontent in Supplementary Section 1 as follows:\nTo rigorously evaluate the contribution of each module in ChrisCell to overall\nperformance, we conducted comprehensive ablation studies. Six model variants were\ncompared: the baseline ChrisCell (500M parameters), ChrisCell with 50M and 100M\nparameters, and three ablated versions—excluding the Gene Ontology (GO)\nembeddings module, the SCD module, and the ZINB loss, respectively.\n`\nFig.1 a–c presents the performance comparison across three clustering tasks, cell\ntype annotation, and gene perturbation prediction. The baseline ChrisCell model\nachieved the highest performance on all datasets and metrics, underscoring the\ncomplementary contributions of each component. The ChrisCell 100M model\noutperformed the 50M variant, yet lagged behind the 500M baseline, indicating that\nscaling model parameters enhances performance. Notably, the ChrisCell 100M model\nalso exceeded the performance of scFoundation—a comparable 100M-parameter\nmodel—across all benchmarks. For instance, ChrisCell 100M attained accuracies of\n0.89, 0.47, and 0.66 on the Pancreas, Zheng, and Segerstolpe datasets, compared to\nscFoundation’s 0.81, 0.43, and 0.46. Similarly, it achieved F1 scores of 0.79, 0.82, and\n0.72, surpassing scFoundation’s 0.68, 0.74, and 0.59. On the Adamson, Norman, and\n`\nDixit perturbation datasets, ChrisCell 100M yielded PCC Delta of 0.612, 0.422, and\n0.399, outperforming scFoundation’s 0.564, 0.388, and 0.256.\nRemoving the ZINB loss resulted in pronounced performance degradation,\nparticularly in clustering and cell type annotation, suggesting increased susceptibility\nto data sparsity and impaired biological representation learning. In contrast, omitting\nthe GO embeddings module had the mildest effect on most tasks, implying a relatively\nmodest role among the components. However, its absence led to more noticeable\ndeclines in gene perturbation prediction, where it attained PCC Delta values of 0.601,\n0.432, and 0.412—comparable to or slightly below the ZINB-ablated model (0.610,\n0.431, 0.409). This indicates that incorporating gene ontology priors enhances the\nmodel’s capacity to capture gene–gene interactions, thereby improving perturbation\nresponse prediction.\nNotably, the model without the SCD module consistently performed worst across\nall settings. Paradoxically, as illustrated in Fig.1d, this variant exhibited the fastest\ntraining convergence. This divergence suggests that while the SCD-ablated model\nconverges rapidly, it fails to capture heterogeneous biological patterns effectively. The\nSCD module thus appears critical for mitigating data heterogeneity and enhancing\ngeneralization, as further discussed in the Section.\n2. Clustering and cell annotation: while the authors show superior performance on\nthe investigated datasets, some more detailed benchmarking (including other baseline\nmodels and single cell datasets) would be beneficial. Recent benchmarking efforts, like\nthe ones from CZI (https://virtualcellmodels.cziscience.com/benchmark-task/cell-\ntype-classification) can provide additional datasets and models for this.\nResponse: Thank you for highlighting this important issue! We conduct comprehensive\ncomparison between ChrisCell with broader baselines on CZI in Supplementary\nSection 2 of the revised manuscript as follows:\nTo comprehensively assess the single-cell representation capacity of ChrisCell, we\nperformed extensive comparisons against eight established methods using the Tabula\nSapiens[1-2] 2.0 (TS2) benchmark under the cz-benchmark framework. Tabula Sapiens\nserves as a reference human cell atlas comprising single-cell transcriptomic data from\n`\nmore than 1.1 million cells across 28 tissues, obtained from 24 healthy donors. The\nupdated Tabula Sapiens 2.0 dataset expands upon the original version, incorporating\nover 500,000 high-quality cells from 27 tissues, with samples derived from both male\n(n = 2) and female (n = 7) individuals. In accordance with cz-benchmark protocols, we\nselected five tissue-specific datasets: Blood, Bone Marrow, Lung, Mammary, and\nThymus. The compared methods include raw gene expression, TranscriptFormer[3],\nscBERT[4], scVI[5], AIDO.Cell[6], UCE[7], scGPT, and GeneFormer[8].\nFollowing the cz-benchmark setup, we conducted cell clustering tasks on the TS2\ndatasets. Cell embeddings were extracted from each method without fine-tuning. A\nnearest-neighbor graph was constructed in the embedding space, and clustering was\nperformed using the Leiden algorithm. The resulting clusters were evaluated against\nground-truth cell type labels using Normalized Mutual Information (NMI),\nAdjusted\nRand Index (ARI), and Silhouette Score (SIL).\nAs illustrated in Fig.2, ChrisCell consistently outperformed all other methods\nacross all datasets, achieving average scores of 0.79 (NMI), 0.50 (ARI), and 0.68 (SIL).\nNotably, raw gene expression exhibited competitive or even superior performance\ncompared to several foundation models. For instance, in the Blood dataset, raw gene\nexpression attained NMI, ARI, and SIL values of 0.66, 0.35, and 0.68, respectively—\nexceeding all baselines except ChrisCell. This suggests that the generalizability of\nexisting foundation models remains constrained despite pretraining on large-scale data,\nhighlighting the challenge posed by data heterogeneity. In contrast, ChrisCell\ndemonstrates robust performance across tissues, indicating its enhanced capability to\nmitigate heterogeneity and improve generalizability.\n`\nWe further evaluated all methods on cell type classification under the same cz-\nbenchmark settings. Embeddings derived from each model were used to train three\nclassifiers: Logistic Regression, K-Nearest Neighbors (KNN), and Random Forest. A\n5-fold cross-validation strategy was applied over five random seeds to ensure statistical\nreliability. Performance was measured using Macro-F1, Precision, and Recall, averaged\nacross folds and classifiers.\nResults in Fig. 3 show that ChrisCell achieves advanced performance across all\ndatasets and metrics, with mean Macro-F1, Precision, and Recall of 0.88, 0.91, and 0.88,\nrespectively. In contrast, raw gene expression yielded substantially lower performance,\nunderscoring how its inherent sparsity limits the ability of conventional classifiers to\ndiscern discriminative features. Foundation models mitigate this issue by condensing\nsparse gene expressions into dense, informative embeddings, thereby improving\nclassification accuracy. The superior performance of ChrisCell underscores its\neffectiveness in overcoming data sparsity and capturing biologically meaningful\nrepresentations conducive to accurate cell type annotation.\n`\nWe further evaluated all methods on cell type classification under the same cz-\nbenchmark settings. Embeddings derived from each model were used to train three\nclassifiers: Logistic Regression, K-Nearest Neighbors (KNN), and Random Forest. A\n5-fold cross-validation strategy was applied over five random seeds to ensure statistical\nreliability. Performance was measured using Macro-F1, Precision, and Recall, averaged\nacross folds and classifiers.\nResults in Fig. 3 show that ChrisCell achieves advanced performance across all\ndatasets and metrics, with mean Macro-F1, Precision, and Recall of 0.88, 0.91, and 0.88,\nrespectively. In contrast, raw gene expression yielded substantially lower performance,\nunderscoring how its inherent sparsity limits the ability of conventional classifiers to\ndiscern discriminative features. Foundation models mitigate this issue by condensing\nsparse gene expressions into dense, informative embeddings, thereby improving\nclassification accuracy. The superior performance of ChrisCell underscores its\neffectiveness in overcoming data sparsity and capturing biologically meaningful\nrepresentations conducive to accurate cell type annotation.\n3. Drug sensitivity prediction: Importantly, the description and the citation to the used\ndrug sensitivity datasets is missing. It is also not clear, which splits the authors used\nin the drug response classification problem. As SCAD performs very badly on this\ntask, some additional baseline model would be needed. For the drug response\nprediction (IC50 regression) task the authors used random split of the train and test\ndata. Random split, especially with IC50 data (Ovchinnikova et al. npj Precision\nOncology 2024) can lead to irrealistic high performance metrics, and can mask real\ndifferences between models. To show the superior performance of ChrisCell, some\nmore realistic train / test splits (new cell lines, new drugs etc.) would be required.\nResponse: Thank you for raising this important point. We have now included the\nrelevant citation (Section 2.2.3) and a detailed description of the drug sensitivity dataset\n(Supplementary Section 10.2). Besides, we include more baselines, including the\nLinear model and Random Forest algorithms.\n`\nRegarding the dataset splitting strategy, we did not employ a random split for the\ndrug sensitivity data. Instead, we adopted a five-fold cross-validation approach: the\nentire dataset was divided into five equally sized folds. For each of the five iterations,\nfour folds were used for training and the remaining one for validation. The final\nperformance metric was obtained by averaging the results across all five folds. This\nmethod is widely accepted in the machine learning community and helps provide a\nrobust and reliable evaluation of model performance. We believe it offers a convincing\nassessment of ChrisCell’s capability in predicting drug sensitivity.\n4. Perturb-seq prediction: Recently several papers showed that simple baselines (like\nlinear or random forest models, some examples in Ahlmann-Eltze et al. Nature\nMethods 2025), and even the simplest train mean model outperform foundation\nmodels for post-treatment RNAseq prediction task. The authors should also compare\ntheir results to these simple models.\nResponse: Thank you for highlighting this important issue! We supplemented the some\nnew baselines, including the linear, GeneCompass and CellFM model, for comparison\nof the Perturb-seq prediction in Section 2.2.4 of the revised manuscript.\n`\nIn Fig.4, we observe that the linear model performs comparably to some deep\nlearning models. Notably, on the Dixit dataset, the MSE, MSE DE, PCC, and PCC DE\nare 0.063, 0.174, 0.82, and 0.74, respectively. These results are comparable to, or even\nbetter than, those of GEARS (0.145, 0.315, 0.68, 0.53), GeneCompass (0.084, 0.182,\n0.7, 0.63), CellFM (0.07, 0.14, 0.85, 0.74), scGPT (0.076, 0.11, 0.73, 0.70), and\nscFoundation (0.066, 0.15, 0.83, 0.74). However, on the Adamson and Norman datasets,\n`\nthe linear model achieves performance that is only comparable to or lower than that of\nGEARS, where these two datasets contain more gene perturbation data. This indicates\nthat while the linear model performs well on smaller datasets, foundation models tend\nto generalize better on larger-scale datasets.\n5. GRN analysis: The identified interactions (like PRSS1 and CTRB1) are probably\nnot gene regulatory, but simple gene co-expression interactions, as none of the\nidentified genes are transcription factors. The authors should compare their gene\ninteraction network with some gene regulatory database, like CollecTRI (Muller-Dott\net al. 2023 Nucleic Acids Research).\nResponse: Thank you for highlighting this important issue! We added the further\nexperiments about the GRN comparison with CollecTRI in Supplementary Section 8,\nas follows:\nTo further validate the capability of ChrisCell-Graph in establishing gene\nregulatory networks and capturing gene–gene interactions, we applied the model to the\nsampled test dataset to infer a unified and generalizable GRN. We focused specifically\non transcription factor (TF) genes, which play pivotal roles in the regulation of key\n`\ndownstream targets. The GRN derived from ChrisCell-Graph—based on attention\nweights—was compared against the curated CollecTRI network for benchmark\nvalidation. We selected two master transcriptional regulators of human skin cells,\nLEF1\nand GATA6, for in-depth analysis. LEF1 (Lymphoid Enhancer-binding Factor 1) acts\nas a central effector of the Wnt/β-catenin signaling pathway. Upon activation, it\ntranslocates to the nucleus where it regulates the expression of target genes critically\ninvolved in cell fate determination, proliferation, and differentiation. GATA6 (GATA\nBinding Protein 6) is a key transcription factor essential for cell lineage specification\nand endoderm-derived organ development, exerting widespread influence over\ndevelopmental and homeostatic processes.Both LEF1 and GATA6 serve as\nfundamental regulators in skin biology, making them ideal candidates for evaluating\nthe biological relevance and accuracy of the GRN inferred by ChrisCell-Graph.\nFig. 7 a and c compare the gene regulatory networks (GRNs) inferred by ChrisCell\nand CollecTRI centered on the transcription factor LEF1. Both models identify several\nestablished downstream targets of the LEF1/β-catenin pathway, including FOXJ1, IL4,\nand DCT, which are known to be directly or indirectly regulated by LEF1. In addition,\nChrisCell-Graph predicts several regulatory interactions not present in CollecTRI, such\nas with GIP, FST, and IDO1—genes previously implicated as functional targets of\nLEF1-mediated Wnt/β-catenin signaling, linking this pathway to metabolic,\nimmunomodulatory, and developmental processes.\nSimilarly, Fig.7b and d show GRNs for GATA6. ChrisCell recovers known\ninteractions documented in CollecTRI, such as with FASLG, THY1, and DUSP16,\nwhich are involved in apoptosis, cell adhesion, and signal transduction and have been\nexperimentally supported as downstream targets of GATA6. Furthermore, ChrisCell-\nGraph identifies novel regulatory candidates—including WNT6, KRAS, and DKK1—\nconsistent with established biological roles: GATA6 acts as a downstream effector of\nboth WNT6 and KRAS signaling while also feedback-regulating the pathway via\nDKK1 induction, thereby fine-tuning processes such as proliferation, differentiation,\nand homeostasis.\nTogether, these results demonstrate that ChrisCell-Graph effectively reconstructs\n`\ngene regulatory networks and uncovers novel, biologically plausible interactions\nbeyond those present in curated databases.\nMinor\n1. What is meant by \"setting the SCD cell embedding at the ”Count” token for the cell\nprompt\". The cell embedding seems to be \"just\" an embedding concatenated to the\ndecoder context and not used in a distinguished way.\nResponse: Thank you for raising this confusion! As we mentioned in Section 2.1, the\nSCD cell embedding is used to represent the global information of the gene expressions\nin the encoder part, where it subsequently concatenated to the decoder context as a\nadditional prompt of global information for predicting the unknown gene expression\nvalues. Finally, it is also used to predict the read depth of this single cell.\n2. Some discussion would be desired on how the \"mainstream\" cell foundation models\nare different from ChrisCell e.g. decoder-only, masked gene expression / gene\npretraining task. It would be fruitful to understand the authors’ intuitions.\nResponse: Thank you for this helpful suggestion! We compare three paradigms of\nsingle-cell foundation models as follows:\nFig.13A, Masked Gene Single-Cell Foundation Models (scBERT, scFoundation):\nMasked single-cell foundation models employ either an encoder-only (scBERT) or\n`\nencoder-decoder (scFoundation) architecture. These models are trained using a masked\ngene modeling objective, wherein a subset of gene expression values is randomly\nmasked. The model receives both masked and unmasked genes as input and is trained\nto reconstruct the masked expression values based on the context provided by the\nunmasked genes. This approach facilitates the learning of gene co-expression patterns\nand enables effective representation learning at the single-cell level.\nFig.13B, Generative Single-Cell Foundation Models (scGPT):Generative single-cell\nfoundation models adopt a decoder-only architecture and are trained with an auto-\nregressive objective. In this framework, the model processes all genes simultaneously\nbut employs causal attention masking to enforce an autoregressive factorization—each\ngene can only attend to preceding genes, thereby sequentially predicting the entire gene\nexpression profile. This design ensures that the prediction for each gene is conditioned\nsolely on the expression values of genes that come before it in the sequence.\nFig.13C, Cell-Gene-Prompt-Based Single-Cell Foundation Models (ChrisCell):\nCell-gene-prompt-based models, such as ChrisCell, utilize an encoder–SCD–decoder\narchitecture and are trained via a cell-gene-prompt (CGP) based zero-inflated negative\nbinomial (ZINB) prediction task. Initially, the encoder module encodes down-sampled\nnon-zero genes into gene embeddings, which are then aggregated via a pooling\noperation to derive a cell-level embedding. These representations are further refined\nthrough a structured contextual decoding (SCD) module. For decoding, all expression\nvalues are initialized to zero. The refined SCD gene embeddings are inserted as prompts\nat corresponding gene positions, and the SCD cell embedding serves as a global prompt.\nThe decoder subsequently predicts the ZINB distribution parameters for all non-\nsampled gene expression values based on these prompted embeddings.\nWe have added this content in Supplementary Section 12 in the revised manuscript.\n3. ChrisCell embeddings were used as a replacement for GO graph of co-expression\ngraph of GEARS?\nResponse: Thank you for this confusion! We still use the GO to build the co-expression\ngraph of GEARS, but we use the ChrisCell embeddings to replace the GEARS’s\nembeddings for gene representation initiation for enhancement .\n`\n4. Based on my knowledge, the Adamson perturb-seq dataset contains only 100\ndifferent perturbations. In Figure 4E, on the left panel there are many more data\npoints. What does each data point represent? A pseudobulked perturbation (as in the\noriginal GEARS evaluation), or individual single cells?\nResponse: Thank you for arising this confusion! A data point in Fig. 4e represent a\nsingle cell in our previous manuscript. In GEARS evaluation, each perturbed cell is\nrandomly assigned with a control cell (unperturbed), where the model predict the\nperturbed gene expression based on the control cell and calculate the difference between\nthe prediction and perturbed cell. For more comprehensively visualizing these\ncomparison, we provide the data point for each perturbation in the Fig.4 of the revised\nmanuscript, where we average difference between the prediction and the perturbed gene\nexpression in the Fig. 4e in the revised manuscript.\n5. The authors write: “Fig. 5b presents the overall performance of ChrisCell and other\nbaseline methods across three perturb-seq datasets, evaluated using four metrics.” -\nin the drug (not perturb-seq) perturbation prediction section.\nResponse: Thank you for pointing this issue! We have corrected this typo in our revised\nmanuscript.\nReference\n[1] Consortium*, T.T.S., Jones, R.C., Karkanias, J., Krasnow, M.A., Pisco, A.O., Quake, S.R.,\nSalzman, J., Yosef, N., Bulthaup, B., Brown, P., et al.: The tabula sapiens: A multiple-organ, single-\ncell transcriptomic atlas of humans. Science 376(6594), 4896 (2022)\n[2] Quake, S.R., Consortium, T.S.: Tabula sapiens reveals transcription factor expression,\nsenescence effects, and sex-specific features in cell types from 28 human organs and tissues. bioRxiv,\n2024–12 (2024)\n`\n[3] Pearce, J.D., Simmonds, S.E., Mahmoudabadi, G., Krishnan, L., Palla, G., Istrate, A.-M.,\nTarashansky, A., Nelson, B., Valenzuela, O., Li, D., et al.: A cross-species generative cell atlas across\n1.5 billion years of evolution: The transcriptformer single-cell model. bioRxiv, 2025–04 (2025)\n[4] Yang, F., Wang, W., Wang, F., Fang, Y., Tang, D., Huang, J., Lu, H., Yao, J.: scbert as a large-\nscale pretrained deep language model for cell type annotation of single-cell rna-seq data. Nature\nMachine Intelligence 4(10), 852–866 (2022)\n[5] Lopez, R., Regier, J., Cole, M.B., Jordan, M.I., Yosef, N.: Deep generative modeling for single-\ncell transcriptomics. Nature methods 15(12), 1053–1058 (2018)\n[6] Ho, N., Ellington, C.N., Hou, J., Addagudi, S., Mo, S., Tao, T., Li, D., Zhuang, Y., Wang, H.,\nCheng, X., et al.: Scaling dense representations for single cell with transcriptome-scale context.\nbioRxiv, 2024–11 (2024)\n[7] Rosen, Y., Roohani, Y., Agarwal, A., Samotorˇcan, L., Consortium, T.S., Quake, S.R., Leskovec,\nJ.: Universal cell embeddings: A foundation model for cell biology. bioRxiv, 2023–11 (2023)\n[8] Cui, Z., Xu, T., Wang, J., Liao, Y., Wang, Y.: Geneformer: Learned gene compression using\ntransformer-based context modeling. In: ICASSP 2024-2024 IEEE International Conference on\nAcoustics, Speech and Signal Processing (ICASSP), pp. 8035–8039 (2024).\nReviewer #3 (Remarks to the Author):\nWang et al. present ChrisCell as a novel single-cell foundation model by combining\nVQ-VAE. The model was shown to outperform other large language models on several\nsingle-cell data analysis tasks. However, the comparative analysis is not\ncomprehensive enough, as only a limited number of baseline methods are included.\nThe introduction of ChrisCell-Graph as an interpreter is interesting, but its actual\ncontribution to the overall interpretability of the model remains unclear. Many\ntechnical details of ChrisCell-Graph are not sufficiently explained and current results\ncan’t be repeated or well evaluated.\nI have the following major concerns:\n1. Lack of Comprehensive Comparison with Existing Methods. The comparative\nstudy with scGPT and scFoundation is narrow. The author need include recent large\n`\nlanguage models, such as UCE[1], GeneCompass[2], and CellFM[3] for comparison.\nResponse: Thank you for raising this helpful suggestion! We conduct more\ncomprehensive comparison with broader methods, including, scBERT, UCE,\nGeneCompass, CellFM, on three single-cell clustering, cell type annotation and gene\nperturbation prediction tasks in Section 2.2.1, 2.2.2 and 2.2.4 of the revised manuscript.\nUpon reviewing these results, we observe that both UCE and CellFM achieve\n`\nperformance that is either comparable to or slightly superior than that of scGPT and\nscFoundation. For instance, the NMI values of UCE across the three datasets are 0.88,\n0.44, and 0.58, respectively, while those of CellFM are 0.89, 0.46, and 0.68. In\ncomparison, scGPT and scFoundation yield NMI values of 0.88, 0.43, and 0.41, and\n0.88, 0.44, and 0.64, respectively. Nevertheless, all these baseline methods still\nunderperform relative to ChrisCell on both cell clustering, cell type annotation and\nproperty prediction tasks. Additionally, we note that GeneCompass and scBERT exhibit\nmoderately lower performance compared to the other baselines, though they still\nsurpass the results obtained using raw gene expression data.\n`\nIn Fig. 4, We have also included several methods for comparison, namely Linear\nmodel, GeneCompass, and CellFM. Among these, GeneCompass demonstrates\nperformance that is either comparable to or slightly inferior to that of scGPT; while\nCellFM achieves results on par with those of scFoundation. Notably, we observe that\nthe linear model performs comparably to some deep learning models. Notably, on the\nDixit dataset, the MSE, MSE DE, PCC, and PCC DE are 0.063, 0.174, 0.82, and 0.74,\nrespectively. These results are comparable to, or even better than, those of GEARS\n(0.145, 0.315, 0.68, 0.53), GeneCompass (0.084, 0.182, 0.7, 0.63), CellFM (0.07, 0.14,\n0.85, 0.74), scGPT (0.076, 0.11, 0.73, 0.70), and scFoundation (0.066, 0.15, 0.83, 0.74).\nHowever, on the Adamson and Norman datasets, the linear model achieves performance\nthat is only comparable to or lower than that of GEARS, where these two datasets\ncontain more gene perturbation data. This indicates that while the linear model performs\nwell on smaller datasets, foundation models tend to generalize better on larger-scale\ndatasets.\n2. Insufficient Justification of Innovation. The core innovation is stated to be the\nintegration of SCD. However, the improvements brought by SCD are unclear, and it\nis also not well demonstrated how SCD addresses the challenge of heterogeneous\ndatasets. The authors should provide a comparison between the encoder-decoder and\nencoder-SCD-decoder architectures under similar parameter settings.\nResponse: Thank you for highlighting this important issue! We conduct comprehensive\nablation studies in Supplementary Section 1 for SCD module, the comparison between\nscFoundation (encoder-decoder 100M) and Chriscell 100M and other modules.\nFurthermore, we extends through discussion for how SCD module overcome the data\nheterogeneity in Supplementary Section 9, as follows:\nAccording to the Fig.1, compared to the strong baseline scFoundation, which also\nconsists of 100M parameters, the ChrisCell 100M model still achieves better\nperformance over scFoundation. For instance, the accuracy of ChrisCell 100M on three\ndatasets (Pancreas, Zheng and Segerstolpe) is 0.89, 0.47, and 0.66, which is obviously\nhigher than scFoundation's 0.81, 0.43, and 0.46; the f1 score of ChrisCell 100M on\nthree datasets is 0.79, 0.82 and 0.72, which also exceeds scFoundation's 0.68, 0.74, and\n`\n0.59; the PCC Delta of ChrisCell 100M on Adamson, Norman and Dixit datasets is\n0.612, 0.422 and 0.399, while the scFoundation is 0.564, 0.388 and 0.256.\nNotably, we find that the model without the SCD module obviously poses the worst\nperformance on all datasets and metrics. However, Fig. 1d illustrates that the model\nwithout the SCD module poses the best training converge among all models. This\nphenomenon provides strong evidence that the model without the SCD module suffers\nfrom data heterogeneity, and the SCD module greatly overcomes this issue.\n`\nThe Supplementary section 9 addresses three fundamental questions:\n(1) What is the origin of data heterogeneity in single-cell data?\n(2) How does data heterogeneity impair the performance of single-cell foundation\nmodels?\n(3) How does the SCD module mitigate data heterogeneity?\nData heterogeneity is a well-recognized challenge across numerous database fields.\nIn single-cell genomics, it arises predominantly from batch effects—technical\nvariations introduced when data are generated across different sources, platforms, or\nexperimental protocols. These discrepancies in processing methods and measurement\nscales result in systematic biases that obscure biological signals.\nWhy is data heterogeneity problematic? When single-cell foundation models are\npre-trained on large-scale datasets comprising multiple batches, they may learn to\nexploit batch-specific artifacts rather than biologically meaningful patterns. This leads\nto information leakage: models can achieve high performance by recognizing batch\n`\nidentifiers, especially since cells within a batch often share cell types or expression\nprofiles. However, such models generalize poorly to unseen data containing new\nbatches, as they fail to adapt to technical variations absent during training.\nTo validate this, we compared the batch correction capability of ChrisCell and an\nSCD-ablated variant using both a training subset (10,000 cells from 452 donors) and a\nheld-out test set. As shown in Fig.5 c–d, the model without SCD achieved an AvgBatch\nscore of 18.33 on the training data, but only 9.86 on the test set—indicating significant\noverfitting to batch effects. This conclusion is further supported by the ablation results\nin Fig.1, where the same model exhibited rapid training convergence yet delivered the\npoorest downstream performance across all tasks.\nFinally, how does the SCD module counteract data heterogeneity? Parallel insights\ncan be drawn from protein language modeling. For instance, AlphaFold2’s predictions\ncan be influenced by template-induced bias, leading to structural inaccuracies. Methods\nsuch as Saprot and ESM3 introduced vector quantization (VQ) to project amino acid\nembeddings into a discrete, unified space, thereby reducing superficial biases while\npreserving essential biological information.\nSimilarly, the SCD module in ChrisCell projects cell and gene embeddings into a\nuniversal, discrete vector space. This process emphasizes principal biological\ncomponents while filtering out batch-specific variations through discretization and\nunification. As a result, the SCD embeddings capture fundamental cellular features\nwithout batch effects. Empirical evidence supports this mechanism: ChrisCell\ndemonstrates consistent AvgBatch scores across both training and test sets (Fig.5c-d),\nand maintains robust performance in ablation studies (Fig.1), confirming the efficacy\nof SCD in overcoming data heterogeneity.\n3. Unclear Description and Application of ChrisCell-Graph. The manuscript provides\na confusing analysis of ChrisCell and ChrisCell-Graph, making it difficult to\nunderstand the relationship and distinctions between the two modules. It is unclear\nwhich tasks use ChrisCell and which use ChrisCell-Graph. For example, Section 2.3.1\nintroduces ChrisCell-Graph, but the identification of marker genes only uses the\nembeddings from ChrisCell; Section 2.3.3 is titled \"Gene Regulatory Network\n`\nAnalysis with ChrisCell-Graph,\" but the content does not mention ChrisCell-Graph\nat all or explain how it is applied in this context.\nResponse: Thank you for highlighting this important issue! We sincerely apologize for\nthis misunderstanding. In Section 2.3.1, we solely utilize ChrisCell's embeddings and\ncell codes to perform marker gene discovery, without incorporating ChrisCell-Graph.\nIn Section 2.3.3, we primarily use ChrisCell-Graph to generate the gene regulatory\nnetwork, while the cell codes produced by ChrisCell are employed to define the cell\nstates. To eliminate any confusion, we have moved the first paragraph of Section 2.3.1\nto Section 2.3 and revised the title of Section 2.3.1 to \"Marker Gene Discovery and\nAnalysis with ChrisCell.\" Furthermore, we have changed the title of Section 2.3.3 to\n\"Gene Regulatory Network Analysis with ChrisCell-Graph\" and fixed some typos,\nclearly indicating ChrisCell-Graph is used in Section 2.3.3.\nBesides, we want to further clarify the relationship and distinctions of ChrisCell\nand ChrisCell-Graph:\nChrisCell is a pretraining model for learning robust cell and gene embeddings, and\ngenerate the “cell codes” by the SCD module. Subsequently, ChrisCell-Graph will\nreceive the cell codes, gene and cell embeddings from ChrisCell to generate graph\nknowledge network for the analysis of downstream tasks.\n4. Please specify the tasks involving fine-tuning and those based solely on embedding\nextraction. Are the same settings used for the comparison with other foundation\nmodels? The comparisons may not be fair if your model used the finetuning while\nother models used only the embedding.\nResponse: Thank you for pointing this issue! In our experiment setting, all prediction\ntasks were performed under the zero-shot setting without fine-tuning. In all analysis\n`\ntasks, we only utilize the embeddings provided by ChrisCell pre-training model(In\naddition to multi-omics tasks, we finetuned ChrisCell model) and train ChrisCell-Graph\non the corresponding dataset. We have added this clarification in the Section 2.1 in the\nrevised manuscript.\n5. The author may consider including more baseline methods of cell type prediction\n(See references https://genomebiology.biomedcentral.com/articles/10.1186/s13059-\n021-02480-2 and https://genomebiology.biomedcentral.com/articles/10.1186/s13059-\n019-1795-z). In addition, scBERT has shown excellent performance in cell type\nannotation. Does ChrisCell also achieve superior or at least comparable performance\nto scBERT in this task?\nResponse: Thank you for providing this good suggession! We conduct comprehensive\ncomparison between ChrisCell with broader baselines on CZI in Supplementary\nSection 2 of the revised manuscript as follows:\nTo comprehensively assess the single-cell representation capacity of ChrisCell, we\nperformed extensive comparisons against eight established methods using the Tabula\nSapiens[1-2] 2.0 (TS2) benchmark under the cz-benchmark framework. Tabula Sapiens\nserves as a reference human cell atlas comprising single-cell transcriptomic data from\nmore than 1.1 million cells across 28 tissues, obtained from 24 healthy donors. The\nupdated Tabula Sapiens 2.0 dataset expands upon the original version, incorporating\nover 500,000 high-quality cells from 27 tissues, with samples derived from both male\n(n = 2) and female (n = 7) individuals. In accordance with cz-benchmark protocols, we\nselected five tissue-specific datasets: Blood, Bone Marrow, Lung, Mammary, and\nThymus. The compared methods include raw gene expression, TranscriptFormer[3],\nscBERT[4], scVI[5], AIDO.Cell[6], UCE[7], scGPT, and GeneFormer[8].\n`\nFollowing the cz-benchmark setup, we conducted cell clustering tasks on the TS2\ndatasets. Cell embeddings were extracted from each method without fine-tuning. A\nnearest-neighbor graph was constructed in the embedding space, and clustering was\nperformed using the Leiden algorithm. The resulting clusters were evaluated against\nground-truth cell type labels using Normalized Mutual Information (NMI), Adjusted\nRand Index (ARI), and Silhouette Score (SIL).\nAs illustrated in Fig.2, ChrisCell consistently outperformed all other methods\nacross all datasets, achieving average scores of 0.79 (NMI), 0.50 (ARI), and 0.68 (SIL).\nNotably, raw gene expression exhibited competitive or even superior performance\ncompared to several foundation models. For instance, in the Blood dataset, raw gene\n`\nexpression attained NMI, ARI, and SIL values of 0.66, 0.35, and 0.68, respectively—\nexceeding all baselines except ChrisCell. This suggests that the generalizability of\nexisting foundation models remains constrained despite pretraining on large-scale data,\nhighlighting the challenge posed by data heterogeneity. In contrast, ChrisCell\ndemonstrates robust performance across tissues, indicating its enhanced capability to\nmitigate heterogeneity and improve generalizability.\nWe further evaluated all methods on cell type classification under the same cz-\nbenchmark settings. Embeddings derived from each model were used to train three\nclassifiers: Logistic Regression, K-Nearest Neighbors (KNN), and Random Forest. A\n5-fold cross-validation strategy was applied over five random seeds to ensure statistical\n`\nreliability. Performance was measured using Macro-F1, Precision, and Recall, averaged\nacross folds and classifiers.\nResults in Fig. 3 show that ChrisCell achieves advanced performance across all\ndatasets and metrics, with mean Macro-F1, Precision, and Recall of 0.88, 0.91, and 0.88,\nrespectively. In contrast, raw gene expression yielded substantially lower performance,\nunderscoring how its inherent sparsity limits the ability of conventional classifiers to\ndiscern discriminative features. Foundation models mitigate this issue by condensing\nsparse gene expressions into dense, informative embeddings, thereby improving\nclassification accuracy. The superior performance of ChrisCell underscores its\neffectiveness in overcoming data sparsity and capturing biologically meaningful\nrepresentations conducive to accurate cell type annotation.\n6. 4.3. In the drug response classification task, SCAD is trained on labeled bulk RNA-\nseq data and unlabeled single cell data. It is not clear how ChrisCell is trained and\nevaluated in this setting.\nResponse: Thank you for pointing this confusion! In fact, instead of training ChrisCell,\nwe leverages ChrisCell to obtain the embeddings of single-cell and bulked data to train\nSCAD models for predicting a cell's sensitivity to the drug. We compare ChrisCell\nperformance with that of scGPT and scFoundation under the same setting. We have\nsupplemented these settings in Section 2.2.3 in the revised manuscript.\n7. In the drug perturbation prediction task, please clarify how the data are split for\ntraining and testing. Additionally, it would be beneficial to include comparisons with\nseveral non-foundational models, such as scVIDR[4] or CellOT[5].\nResponse: Thank you for raising this important point. We added these baselines,\nincluding scVIDR and CellOT for more comprehensive comparison in the Section 2.25\nin the revised manuscript.\n`\nFig. 5b presents the overall performance of ChrisCell and other baseline methods on\nthe OP3 datasets, evaluated using four metrics. The results indicate that relying solely\non raw gene expression data yields suboptimal performance, with MSE and MSE DE\nvalues exceeding 1.25 and 4, respectively, while the PCC and PCC DE values fall below\n0.66 and 0.46. This suggests that predicting using raw gene expression data does not\nadequately align with the ground truth. Besides, we noticed that the non-foundation\nmodels can achieve comparable or better performance over some foundation models.\n`\nFor instance, the MSE, MSE DE, PCC and PCC DE of CellOT are 0.247, 0.763, 0.901\nand 0.782, which are superior than scGPT(0.887, 2.333, 0.853, 0.678) and\nscFoundation(0.312, 1.163, 0.897, 0.800). However, ChrisCell still achieves better\nperformance over other methods. The representation derived from ChrisCell reduces\nMSE and MSE DE to 0.087 and 0.233, respectively, while improving PCC and PCC\nDE to 0.965 and 0.945. These results demonstrate a significant enhancement in\nprediction quality, with ChrisCell outperforming other robust baselines, such as\nscVIDR and scFoundation, across all metrics.\nFig. 5c compares the Pearson delta and Pearson delta DE of ChrisCell with five\nbaseline methods on the OP3 dataset. The average Pearson delta and Pearson delta DE\nfor ChrisCell are 0.478 and 0.677, respectively. In comparison, raw gene expression\nyields values of 0.12 and 0.223, while scGPT reports 0.235 and 0.356, scVIDR is 0.212,\n0.343, CellOT reports 0.268, 0.391, and scFoundation yields 0.267 and 0.388. ChrisCell\nexhibits superior performance in both Pearson delta metrics, indicating that its gene\nrepresentation captures more co-expression information, thereby facilitating more\naccurate predictions.\n8. The meaning of Supplementary Fig. 3a is unclear. Please clarify what this figure\nrepresents and explain how it supports the claim of validating the effectiveness of\nChrisCell-Graph’s attention weights.\nResponse: Thank you for raising this important point. Indeed, we have provided a\nclarification regarding Supplementary Fig. 3a in Section 2.3.2 of the revised manuscript,\nwhich reads as follows:\n“We calculated the Mean Square Distance between the vectors corresponding to\ndifferent states and visualized this information in the distance confusion matrix\npresented in Supplementary Fig.3a. This visualization supports our earlier conclusions.\nSpecifically, states 3302, 3430, and 3174 exhibit the closest distances among all states,\nwhile states 1474 and 1490 also display relatively close distances to one another.\nNotably, the distance between states 3174 and 3302 is closer than that between either\nof these states and 3430, corroborating the findings illustrated in Fig.7c. We posit that\nthis characterization provided by ChrisCell-Graph offers a more detailed and precise\n`\nquantization method for measuring the distances between cell states.” Unfortunately,\nin our previous manuscript, the incorrect figure was referred. We apologize for this\noversight and have now corrected the reference in Section 2.3.2 of the revised\nmanuscript.\n9. The process for generating Figure 7d is unclear and requires further explanation.\nResponse: Thank you for raising this important point. For generating Figure 7d, we\nobtain the corresponding SCD embeddings of each cell, performing the UMAP\nvisualization. Subsequently, according to the attention weights of ChrisCell-Graph,\nwe\nannotate the state gene(the gene with the highest attention weight) to each cell state in\nthis UMAP visualization. We have supplemented these contents in the Section 2.3.2\n10. The term \"language tokenizer\" is introduced but not sufficiently explained in the\nmanuscript.\nResponse: Thank you for pointing this limitation. The language tokenizer used in our\npaper is introduced by BioGPT. It include 42,383 tokens and can convert the text of\nsingle cell properties description into one or a series of tokens (intergers), effectively\nembedding these properties within a unified graph framework. We have supplemented\nthese content in the Section 2.3.4.\n11. In Figure 9a-d, the meaning of the x-axis is unclear. Additionally, please clarify\nhow the attention score is calculated and what entities are involved in this calculation.\nResponse: Thank you for pointing this limitation. In Figure 9a-d, the x-axis is the\nnumber of the cell state generated by ChrisCell. The attention score is calculated by the\ngraph attention module of the ChrisCell-Graph. Given the cell embedding h and the\ni\nproperty embedding h, the attention module firstly compute the coefficients e ,\nj ij\nrepresenting the importance of property j's features for cell i:\nThese coefficients are normalized using the softmax function to ensure they sum to 1\nfor cell i:\n`\nThe 𝛼 is the attention weights of the property j for cell i, where the 𝑁 represent the\n𝑖𝑗 𝑖\nproperty set of cell i. These calculation details are listed in the Section 4.4.2 of our\nmanuscript.\n12. In section 2.3.5, the procedure for fine-tuning on ATAC-seq data should be\ndescribed in detail.\nResponse: Thank you for your helpful suggestion! We provided the detailed fine-tuning\nof ChrisCell and the training of ChrisCell-Graph on ATAC-seq data in Supplementary\nSection 11.5:\n“We began by fine-tuning the ChrisCell model on scATAC-seq data. Given the\ncharacteristically high dimensionality and sparsity of scATAC-seq datasets, we\nrandomly sampled a maximum of 4096 peaks as input to the encoder to obtain\nStructured Cell Representation (SCD) embeddings. To mitigate the challenges posed by\ndata sparsity, we sampled an equal number of non-accessible regions for zero-inflated\nnegative binomial (ZINB) distribution prediction. The fine-tuning procedure employed\nthe following hyperparameters: a learning rate of 1e-4, learning rate decay steps of\n1000, 1000 training epochs, and a batch size of 8.\nFollowing the fine-tuning of ChrisCell, we proceeded to train the ChrisCell-Graph\nmodel using SCD embeddings derived from the fine-tuned ChrisCell. This phase utilized\na multi-omics dataset—specifically, the mouse cerebral cortex dataset generated by\nSNARE-seq. The model architecture hyperparameters remained consistent with\nprevious configurations, including a hidden embedding dimension of 128 and a single\ngraph attention layer. Training was conducted with a learning rate of 5e-4, learning\nrate decay steps of 1000, 1000 epochs, and a batch size of 128. Optimization was\nperformed using the AdamW optimizer coupled with a OneCycle learning rate\nscheduler.”\nOther issues:\n13. Most existing large models use MSE as the reconstruction loss, while the authors\nuse ZINB loss. What is the specific benefit of ZINB over MSE? It would be better to\nprovide an ablation experiment comparing the two losses to clarify the impact.\nResponse: Thank you for pointing this limitation! We conducted a thorough ablation\n`\nstudy for ZINB loss in Supplementary Section 1 in the revised manuscript, as follows:\nRemoving the ZINB loss resulted in pronounced performance degradation,\nparticularly in clustering and cell type annotation, suggesting increased susceptibility\nto data sparsity and impaired biological representation learning. In contrast, omitting\nthe GO embeddings module had the mildest effect on most tasks, implying a relatively\nmodest role among the components. However, its absence led to more noticeable\ndeclines in gene perturbation prediction, where it attained PCC Delta values of 0.601,\n0.432, and 0.412—comparable to or slightly below the ZINB-ablated model (0.610,\n0.431, 0.409). This indicates that incorporating gene ontology priors enhances the\nmodel’s capacity to capture gene–gene interactions, thereby improving perturbation\nresponse prediction.\n`\nIn addition, we discuss the contribution of ZINB predictor and loss function in\nSection 4.4 of the revised manuscript: As prior works supported[9,10], ZINB predictor\neffectively addresses data sparsity by combining a zero-inflation model with a negative\nbinomial distribution. It first models additional zeros through the zero-inflation\ncomponent, enabling the model to capture the high proportion of zeros in real data\nrather than treating them solely as outcomes of the negative binomial distribution. The\nnegative binomial component models the variability and discreteness of non-zero count\ndata, providing a more accurate description of sparse data. This combination makes\nZINB particularly well-suited for handling highly sparse data in fields such as gene\nexpression, thereby enhancing the model's fitting capability and predictive accuracy.\n14. There is an unexpected red wavy line appearing in Supplementary Figure 1.\nResponse: Thank you for pointing this error. We have corrected this issue of the\nSupplementary Figure 1 in the Supplementary Figure 4 of the revised manuscript.\n15. In Section 2.3.2 (last paragraph) and Section 2.3.3 (third paragraph), figure\nreferences are incorrect.\nResponse: Thank you for pointing this error. We have corrected these issue of the\nSection 2.3.2 and 2.3.3 in the revised manuscript.\n16. Some biological insights require literature support, such as the marker genes\nshown in Figure 6c.\nResponse: Thank you for this insightful suggestion. In response, we have incorporated\nadditional citations pertaining to marker genes in the revised manuscript, as provided\nin references[47-64] of the revised manuscript.\nReference\n[1] Consortium*, T.T.S., Jones, R.C., Karkanias, J., Krasnow, M.A., Pisco, A.O., Quake,\nS.R., Salzman, J., Yosef, N., Bulthaup, B., Brown, P., et al.: The tabula sapiens: A multiple-\norgan, single-cell transcriptomic atlas of humans. Science 376(6594), 4896 (2022)\n[2] Quake, S.R., Consortium, T.S.: Tabula sapiens reveals transcription factor expression,\nsenescence effects, and sex-specific features in cell types from 28 human organs and tissues.\n`\nbioRxiv, 2024–12 (2024)\n[3] Pearce, J.D., Simmonds, S.E., Mahmoudabadi, G., Krishnan, L., Palla, G., Istrate, A.-\nM., Tarashansky, A., Nelson, B., Valenzuela, O., Li, D., et al.: A cross-species generative\ncell atlas across 1.5 billion years of evolution: The transcriptformer single-cell model.\nbioRxiv, 2025–04 (2025)\n[4] Yang, F., Wang, W., Wang, F., Fang, Y., Tang, D., Huang, J., Lu, H., Yao, J.: scbert as a\nlarge-scale pretrained deep language model for cell type annotation of single-cell rna-seq\ndata. Nature Machine Intelligence 4(10), 852–866 (2022)\n[5] Lopez, R., Regier, J., Cole, M.B., Jordan, M.I., Yosef, N.: Deep generative modeling\nfor single-cell transcriptomics. Nature methods 15(12), 1053–1058 (2018)\n[6] Ho, N., Ellington, C.N., Hou, J., Addagudi, S., Mo, S., Tao, T., Li, D., Zhuang, Y., Wang,\nH., Cheng, X., et al.: Scaling dense representations for single cell with transcriptome-scale\ncontext. bioRxiv, 2024–11 (2024)\n[7] Rosen, Y., Roohani, Y., Agarwal, A., Samotorˇcan, L., Consortium, T.S., Quake, S.R.,\nLeskovec, J.: Universal cell embeddings: A foundation model for cell biology. bioRxiv,\n2023–11 (2023)\n[8] Cui, Z., Xu, T., Wang, J., Liao, Y., Wang, Y.: Geneformer: Learned gene compression\nusing transformer-based context modeling. In: ICASSP 2024-2024 IEEE International\nConference on Acoustics, Speech and Signal Processing (ICASSP), pp. 8035–8039 (2024).\n[9] Kalfon, J., Samaran, J., Peyre, G., Cantini, L.: scprint: pre-training on 50million cells\nallows robust gene network predictions. Nature Communications 16(1), 3607 (2025)\n[10] Lopez, R., Regier, J., Cole, M.B., Jordan, M.I., Yosef, N.: Deep generative modeling\nfor single-cell transcriptomics. Nature methods 15(12), 1053-1058 (2018)\nReviewer #4 (Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed\nreports. This is part of the Nature Communications initiative to facilitate training in\npeer review and to provide appropriate recognition for Early Career Researchers who\nco-review manuscripts.\n`\nReviewer #4 (Remarks on code availability):\nThe code is comprehensive, but it lacks guidance on testing for each downstream task.\nResponse: Thank you for this helpful suggestion! In fact, we have provided all source\ncodes and README files in our repository in each directory. But we did not mention\nthat in our main README, we apologize for this negligence and have added this note\nto our README.\n`"}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have responded to and resolved the majority of concerns raised in the initial review. However, regarding the second point raised in my previous comments, I still require clarification on two key issues: First, the biological content corresponding to the referenced numerical state identifiers (e.g., 3174, 1474) needs explicit definition in the manuscript. Second, the selection criteria for choosing these particular states for detailed presentation rather than others must be clearly justified. Furthermore, it would be valuable to explore whether the analysis could be extended to all states in a comprehensive manner - either through global assessment or grouped comparisons - to provide proper context for the selected examples.", "reviewer_2": "1. Ablation studies: the authors answered our questions, the ablation studies now can help the readers to understand the effect of different architectural choices.\n2. Clustering and cell annotation: The authors answered our questions.\n3. Drug sensitivity prediction: Thank the authors for the clarification. However, we still think that the current splitting / cross-validation setup is not optimal to measure ChrisCell performance. If we understand correctly, in the IC50 prediction task the authors use drug representations and cell representations (coming from ChrisCell and other methods) to predict drug sensitivity. Train and test samples were randomly split (even if the authors used a 5-fold CV setup), so each cell line can be present both in the train and the test set (with different drugs of course). A much better approach to test drug sensitivity prediction would be to leave out entire cell lines from the train test. This would allow to test how the cell representation (coming from ChrisCell) helps generalization.\n4. Perturb-seq prediction: There is some discrepancy between the authors results and the Ahlmann-Eltze paper results. E.g.: In the Ahlmann-Eltze paper, the linear model clearly outperforms scGPT on the Adamson dataset (Extended Data Figure 8), which is not the case in this manuscript. Also a more detailed description of the baseline methods would be beneficial.\n5. GRN analysis: The authors show GRN coming from ChrisCell and Collectri. However these GRN plots are very hard to interpret - some evaluation metrics (accuracy, overlap coefficient, Jaccard index etc.) and comparison to some “random” model would be required. Also, the authors explain in their answer that some of the ChrisCell identified interactions are biologically important - here citations for the source of these informations would be required. Minor: They refer to the ChrisCell - Collectri comparison a Figure 7 in their answer - but it is Figure 9 in Supplementary?", "reviewer_3": "We thank authors for the detailed responses. The manuscript has improved and most of the earlier concerns have been addressed. However, a few key aspects of ChrisCell-Graph’s motivation and implementation remain unclear. Clarifying the points below would further strengthen the paper.\n1. The relationship between ChrisCell and ChrisCell-Graph is confusing. ChrisCell-Graph does not show a clear advantage over ChrisCell. In many analyses, it seems that the same results could be achieved using ChrisCell’s embeddings and attention weights alone.\na) Figures 6b and 7b both show that both models can cluster cells of the same type well, but it leaves an important question: how do the cell embeddings produced by the two methods differ? A direct comparison would help clarify their respective strengths.\nb) Figures 6c–e and Figure 7e involve two different approaches for identifying key genes: max-pooling from ChrisCell and attention from ChrisCell-Graph. What are the differences between these two methods? Under what circumstances should ChrisCell be used, and when is ChrisCell-Graph the better choice?\n2. In section 2.3.3, the authors state “Two genes are defined as connected if they share the same code”. How are genes connected to cell codes in the model? In addition, models such as scGPT construct GRNs directly by computing the similarity between gene embeddings. What advantages does the proposed approach offer compared to scGPT?\n3. The training process for ChrisCell-Graph remains unclear. Please clarify how ChrisCell-Graph is trained, including whether the training objectives are supervised or unsupervised. How the cell graph is constructed? In particular, how are the edges between cells defined? In Section 2.3, the authors state that ‘ChrisCell generates both cell and gene embeddings for ChrisCell-Graph’. Does the graph in ChrisCell-Graph actually include genes as nodes? This point is not evident from Figure 6a.\n4. Since ChrisCell is trained on single-cell transcriptomic (scRNA-seq) data, how can it be directly fine-tuned on scATAC-seq data? The encoding spaces of these two data types are inherently different.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "NCOMMS-25-59714-T response\nWe thank all reviewers for the insightful suggestions and comments to help us improve\nthe manuscript. Below are our point-by-point responses to address all the concerns\nraised by reviewers. All the changes made in the main manuscript are highlighted in\nyellow.\nOverview of the Revisions:\n1. We gave a explicit definition for the cell states provided by ChrisCell\nTo address the comment from Reviewer 1 (Question 1), we have now provided an\nexplicit biological definition for the cell states characterized by ChrisCell, grounded\nin citations from the relevant literature. Furthermore, we clarified the rationale\nbehind the selection of these specific cell states and underscored that the analytical\ncapabilities of ChrisCell and ChrisCell-Graph are generalizable.\n2. We conducted a detailed description for baseline methods.\nIn response to Reviewer 2's Question 4, we provided an explanation for the\ndiscrepancies between our results and those in the Ahlmann-Eltze paper.\nFurthermore, we provided a detailed description for baseline methods in\nSupplementary Section 14.\n3. We gave a metric to assess the comparison of the GRNs constructed by\nChrisCell and ColleTRI and add a random model for comparison.\nIn response to Reviewer 2's Question 5, we have provided the overlap coefficients\nbetween the GRNs constructed by ChrisCell and ColleTRI for assessment.\nFurthermore, we supplemented this analysis with a random model for additional\ncomparison.\n4. We supplemented some further describtion for the construction of graph and\nthe training process of ChrisCell-Graph.\nIn response to Question 2, 3 of Reviewer 3, we have supplemented the describtion\nfor the graph construction and the training process of ChrisCell-Graph to to prevent\nany confusion.\n`\nRespond List to Comments:\nReviewer #1 (Remarks to the Author):\nThe authors have responded to and resolved the majority of concerns raised in the\ninitial review. However, regarding the second point raised in my previous comments,\nI still require clarification on two key issues: First, the biological content\ncorresponding to the referenced numerical state identifiers (e.g., 3174, 1474) needs\nexplicit definition in the manuscript. Second, the selection criteria for choosing these\nparticular states for detailed presentation rather than others must be clearly justified.\nFurthermore, it would be valuable to explore whether the analysis could be extended\nto all states in a comprehensive manner - either through global assessment or grouped\ncomparisons - to provide proper context for the selected examples.\nResponse: Thank you for raising this important point. Regarding the first question, we\nused oligodendrocyte cell states as an example to provide a clear biological definition.\nAs supported by several studies[1,2], oligodendrocytes exhibit considerable\nheterogeneity, comprising distinct subpopulations and dynamic states associated with\ndifferent stages of lineage progression and potentially reflecting diverse functional roles.\nThis concept dates back to Río-Hortega’s early morphological classification of four\noligodendrocyte types[3], which has since been further refined in modern single-cell\nstudies[4,5].\nTo simplify this complexity, recent literature[1,2] categorizes oligodendrocyte\ndevelopment into three main states:\nNewly Differentiated Oligodendrocyte: This state corresponds to cells that have\njust exited the oligodendrocyte precursor cell (OPC) stage. They begin expressing\nclassic markers such as myelin basic protein (MBP) and ANK3, display a simple\nbranched morphology, and have not yet initiated extensive myelination. In our study,\nthis corresponds to state 3302.\nMyelinating-Competent Oligodendrocyte: An intermediate, highly motile state\nwhere cells extend multiple processes to survey and establish contact with axonal\n`\ntargets. This represents a transitional stage prior to active myelination, corresponding\nto state 3174 in our work.\nMature Myelinating Oligodendrocyte: The terminally differentiated state\ncharacterized by the formation of compact, multi-layered myelin sheaths around axons.\nThese cells are metabolically active and essential for long-term axonal integrity and\nrapid signal conduction. This state aligns with state 3430 in our analysis.\nRegarding the second question, we emphasize that diverse cell types—including\noligodendrocytes, central nervous system macrophages, and others in our paper—\nexhibit complex and dynamic functional states that are challenging to delineate[1,2,6,7].\nComputational approaches such as ChrisCell and ChrisCell-Graph provide a\nquantitative and traceable framework for identifying these states in a biologically\nmeaningful way.\nMoreover, the analytical capability of ChrisCell-Graph is universally applicable\nacross all cell states. By constructing knowledge graphs and modeling relationships\namong genes and cells within each state, ChrisCell-Graph detects key state-specific\ngenes. Through functional annotation of these genes or integration with other\nexperimental tools, these computationally derived states become biologically\ninterpretable and validated. We have added this clarification to the Supplementary\nSection 6.\n[1] van Bruggen, David, Eneritz Agirre, and Goncalo Castelo-Branco. \"Single-cell\ntranscriptomic analysis of oligodendrocyte lineage cells.\" Current opinion in neurobiology 47\n(2017): 168-175.\n[2] Huang, Wei, et al. \"Origins and proliferative states of human oligodendrocyte precursor\ncells.\" Cell 182.3 (2020): 594-608.\n[3] H. Lassmann The birth of oligodendrocytes in the anatomical and neuropathological\nliterature: the seminal contribution of Pio Del Rio-hortega, 1921\n[4] E.S. Anderson, C. Bjartmar, C. Hildebrand Myelination of prospective large fibres in chicken\nventral funiculus J Neurocytol, 29 (2000), pp. 755-764\n[5] C. Bjartmar, C. Hildebrand, K. Loinder Morphological heterogeneity of rat oligodendrocytes:\n`\nelectron microscopic studies on serial sections\n[6] Prinz, Marco, et al. \"Microglia and central nervous system–associated macrophages—from\norigin to disease modulation.\" Annual review of immunology 39.1 (2021): 251-277.\n[7] Goldmann, Tobias, et al. \"Origin, fate and dynamics of macrophages at central nervous\nsystem interfaces.\" Nature immunology 17.7 (2016): 797-805.\nReviewer #2 (Remarks to the Author):\n1. Ablation studies: the authors answered our questions, the ablation studies now can\nhelp the readers to understand the effect of different architectural choices.\n2. Clustering and cell annotation: The authors answered our questions.\n3. Drug sensitivity prediction: Thank the authors for the clarification. However, we\nstill think that the current splitting / cross-validation setup is not optimal to measure\nChrisCell performance. If we understand correctly, in the IC50 prediction task the\nauthors use drug representations and cell representations (coming from ChrisCell and\nother methods) to predict drug sensitivity. Train and test samples were randomly split\n(even if the authors used a 5-fold CV setup), so each cell line can be present both in\nthe train and the test set (with different drugs of course). A much better approach to\ntest drug sensitivity prediction would be to leave out entire cell lines from the train\ntest. This would allow to test how the cell representation (coming from ChrisCell)\nhelps generalization.\nResponse: Thank you for raising this important point! If we understand correctly, your\nconcern is that ChrisCell is trained on cell line data and therefore may not adequately\ndemonstrate its generalization capabilities. However, we would like to clarify that, as\nstated in Section 2.1: “All prediction tasks were performed under the zero-shot setting\nwithout fine-tuning.” ChrisCell itself was not trained on any cell line data used in our\nevaluation. We only trained the downstream module of SCAD to predict IC50 values\nbased on the cell representations provided by ChrisCell. Therefore, we believe our\nexperimental settings address your concern and meet the requirements for\ndemonstrating generalization.\n`\n4. Perturb-seq prediction: There is some discrepancy between the authors results and\nthe Ahlmann-Eltze paper results. E.g.: In the Ahlmann-Eltze paper, the linear model\nclearly outperforms scGPT on the Adamson dataset (Extended Data Figure 8), which\nis not the case in this manuscript. Also a more detailed description of the baseline\nmethods would be beneficial.\nResponse: Thank you for highlighting this important discrepancy! We have identified\nkey differences between the experimental settings in the Ahlmann-Eltze paper and our\nstudy. In the Ahlmann-Eltze paper, the approach was: 'We fine-tuned the models on all\n100 single perturbations and on 62 of the double perturbations and assessed the\nprediction error on the remaining 62 double perturbations.' In contrast, our experimental\nsetting extracts gene representations from these models and trains only the downstream\nmodule of GEARS to predict gene expression perturbations.We believe the discrepancy\nbetween the results of the Ahlmann-Eltze paper and ours stems from two key\nperspectives:\nPotential fine-tuning inconsistencies: Foundation models such as scGPT and\nscFoundation did not provide official training code for fine-tuning. Consequently, the\nfine-tuning implementation in the Ahlmann-Eltze paper may have introduced\nmethodological biases that could affect the results. In our approach, we directly extract\ngene representations from these pre-trained models without fine-tuning, thereby\navoiding such potential biases.\nGeneralization capability: Direct fine-tuning of foundation models on downstream\ndatasets, as done in the Ahlmann-Eltze paper, may significantly compromise the\ngeneralization power of these models, potentially leading to suboptimal performance.\nIn contrast, our approach follows the same setting as the original scFoundation paper\nby using only the gene embeddings from foundation models as gene representations.\nThis methodology better preserves the models' generalization capabilities and avoids\noverfitting to specific downstream tasks.\nThis design choice ensures a fairer and more robust evaluation of the inherent\nrepresentation quality of these foundation models.\nRegarding to the description of baseline methods, we add these contents in\n`\nSupplementary Section 14 as follow:\nGEARS is a generative modeling framework specifically designed for predicting\ntranscriptional responses to genetic perturbations. It leverages a graph-based\narchitecture that incorporates prior knowledge of gene interactions, such as pathway\ninformation or protein-protein interaction networks. This allows GEARS to not only\npredict the effects of single-gene perturbations with high accuracy but also to generalize\nto unseen combinatorial perturbations, making it a powerful in-silico tool for\nhypothesizing about genetic interactions and their outcomes on the transcriptome.\nGeneCompass is a foundational model tailored for spatial transcriptomics data,\ntreating spatial gene expression patterns as a language to be learned. By pre-training on\nvast amounts of spatial data, it captures the complex relationships between gene\nexpression and physical location within a tissue. This enables GeneCompass to perform\na variety of zero-shot or few-shot tasks, such as imputing missing genes, deconvoluting\ncell types, and inferring spatial gene-gene interactions, all through prompt-based\ninference without the need for task-specific retraining.\nCellFM positions itself as a universal cell simulator by learning a foundational\nrepresentation from over 50 million single cells. Its core innovation lies in a \"prompt-\nprediction\" framework where users can input a cell's state and a text-based instruction\n(e.g., \"predict the response to drug X\"). This design allows CellFM to perform versatile\nzero-shot and few-shot predictions across diverse tasks and even across species, aiming\nto generalize to novel biological questions that were not explicitly part of its training\ndata.\nscGPT is a foundational model that applies the Transformer architecture, the\nbackbone of modern large language models, to single-cell biology. It conceptualizes\ngenes as \"words\" and a cell's expression profile as a \"sentence,\" then pre-trains on\nmillions of cells using a masked gene prediction objective. This process allows scGPT\nto learn deep, context-aware relationships between genes. The pre-trained model can\nthen be adapted through fine-tuning or prompting for a wide range of downstream tasks,\nincluding batch correction, cell type annotation, and perturbation prediction,\nestablishing a unified paradigm for single-cell analysis.\n`\nscFoundation is a large-scale generative pre-trained model developed for single-\ncell transcriptomics. Similar in spirit to scGPT, it is built on a Transformer architecture\nand trained on a 50 million of single-cell data to reconstruct masked gene expressions.\nIts key strength is its scale, both in terms of model parameters and the diversity of\ntraining data, which empowers it to capture a rich and robust representation of cell states.\nscFoundation serves as a potent base model that can be efficiently fine-tuned for\nspecialized tasks like predicting dose-responsive gene expressions and enhancing the\nanalysis of rare cell types.\n5. GRN analysis: The authors show GRN coming from ChrisCell and Collectri.\nHowever these GRN plots are very hard to interpret - some evaluation metrics\n(accuracy, overlap coefficient, Jaccard index etc.) and comparison to some “random”\nmodel would be required. Also, the authors explain in their answer that some of the\nChrisCell identified interactions are biologically important - here citations for the\nsource of these informations would be required. Minor: They refer to the ChrisCell -\nCollectri comparison a Figure 7 in their answer - but it is Figure 9 in Supplementary?\nResponse: Thank you for raising this important issue. We have calculated the overlap\ncoefficient between the GRNs inferred by ChrisCell and the reference network from\nCollecTRI. To provide context, we also included GRNs generated by a random model\nin Supplementary Fig. 9. Panels a and c of this figure compare the GRNs centered on\nthe transcription factor LEF1, revealing an overlap coefficient of 0.48. Similarly, panels\nb and d show the GRNs for GATA6, with an overlap coefficient of 0.46.\n`\nFurthermore, we have added citations (Supplementary References 46-48, 52-53,\n57-59) to support the biological relevance of the regulatory interactions identified by\nChrisCell-Graph. We sincerely apologize for our oversight in not making this clearer\ninitially.\nReviewer #3 (Remarks to the Author):\nWe thank authors for the detailed responses. The manuscript has improved and most\nof the earlier concerns have been addressed. However, a few key aspects of ChrisCell-\nGraph’s motivation and implementation remain unclear. Clarifying the points below\nwould further strengthen the paper.\n`\n1. The relationship between ChrisCell and ChrisCell-Graph is confusing. ChrisCell-\nGraph does not show a clear advantage over ChrisCell. In many analyses, it seems\nthat the same results could be achieved using ChrisCell’s embeddings and attention\nweights alone.\na) Figures 6b and 7b both show that both models can cluster cells of the same type\nwell, but it leaves an important question: how do the cell embeddings produced by the\ntwo methods differ? A direct comparison would help clarify their respective strengths.\nb) Figures 6c–e and Figure 7e involve two different approaches for identifying key\ngenes: max-pooling from ChrisCell and attention from ChrisCell-Graph. What are\nthe differences between these two methods? Under what circumstances should\nChrisCell be used, and when is ChrisCell-Graph the better choice?\nResponse: Thank you for raising this crucial point. We would like to clarify the\ndistinction between ChrisCell and ChrisCell-Graph as follows:\nChrisCell is a single-cell foundation model pre-trained on massive single-cell data\nto generate informative and high-quality representations of cells and genes, as well as\ncell codes. However, as a large model with extensive parameters, ChrisCell is inherently\ncomplex and challenging to interpret. Additionally, it cannot directly model multi-\nmodal data.\nChrisCell-Graph is a lightweight, plug-and-play module that takes the cell and gene\nrepresentations from ChrisCell and integrates multi-modal data—such as cells, genes,\nchromatin accessibility, and cell properties—as nodes into a shared knowledge graph.\nIt then learns the relationships between these nodes. Thanks to its lightweight and\nflexible design, ChrisCell-Graph can be easily and rapidly applied to custom datasets\nfor biological knowledge discovery and offers explicit interpretability.\nRegarding question (a), the key difference in cell embeddings between the two\nmethods is that ChrisCell embeddings encode only gene expression information,\nwhereas ChrisCell-Graph embeddings incorporate gene expression, cell properties, and\ncell-cell communication.\nFor question (b), the difference lies in how key genes are identified. In Figure 6c–\n`\ne, marker genes for a specific cell type are derived by calculating the frequency of a\ngene’s channel in the max-pooled ChrisCell embeddings, relying solely on gene\nexpression within individual cells. In contrast, the state genes in Figure 7e are\nidentified using ChrisCell-Graph, which integrates not only gene expression from a\nsingle cell but also cell communication signals from neighboring cells. This approach\nreflects that cell states are influenced by both intrinsic gene expression and the external\ncellular environment.\nIn summary, ChrisCell is sufficient for tasks that aim to identify cell-type marker\ngenes without considering cell communication or other properties. However, for the\nidentification of state-specific genes that define a cell state, ChrisCell-Graph is the\nappropriate tool.\n2. In section 2.3.3, the authors state “Two genes are defined as connected if they share\nthe same code”. How are genes connected to cell codes in the model? In addition,\nmodels such as scGPT construct GRNs directly by computing the similarity between\ngene embeddings. What advantages does the proposed approach offer compared to\nscGPT?\nResponse: Thank you for pointing out this confusion! As noted in Section 2.1—“We\ninput the gene embeddings into an SCD module to obtain SCD gene embeddings”—\nthe gene embeddings do indeed receive a code from the SCD module. However, since\nthese codes did not demonstrate significant importance in our training or downstream\ntasks, we did not emphasize them in the original text. We apologize for this omission\nand have added clarification in Section 2.3.3.\nRegarding GRN construction compared to scGPT, our method offers several\nadvantages. scGPT infers GRNs based on gene embedding similarity, which heavily\ndepends on the choice of threshold to determine whether two genes are connected. As\na result, different thresholds lead to different networks, making the GRN uncertain and\nsensitive to this parameter. In contrast, ChrisCell constructs a deterministic GRN by\nlearning shared biological patterns from large-scale pretraining, providing a more stable\nand principled approach.\n3. The training process for ChrisCell-Graph remains unclear. Please clarify how\n`\nChrisCell-Graph is trained, including whether the training objectives are supervised\nor unsupervised. How the cell graph is constructed? In particular, how are the edges\nbetween cells defined? In Section 2.3, the authors state that ‘ChrisCell generates both\ncell and gene embeddings for ChrisCell-Graph’. Does the graph in ChrisCell-Graph\nactually include genes as nodes? This point is not evident from Figure 6a.\nResponse: Thank you for raising this important point. The training of ChrisCell-Graph\nis a self-supervised task. It takes precomputed cell and gene embeddings from ChrisCell,\nintegrates cell property data to construct a knowledge graph, and is trained by\nminimizing the reconstruction loss of the original gene expression. The cell graph is\nbuilt using a k-nearest neighbors (kNN) algorithm based on the lowest gene expression\ndifference between cells, a method supported by prior work such as scGCOT [1].\nFurthermore, the graph in ChrisCell-Graph incorporates genes as nodes. A triangle in\nFigure 6a can represent a gene as a node to connect its corresponding cell. We have\nadded these clarifications in Supplementary Section 12.4.\n4. Since ChrisCell is trained on single-cell transcriptomic (scRNA-seq) data, how can\nit be directly fine-tuned on scATAC-seq data? The encoding spaces of these two data\ntypes are inherently different.\nResponse: Thank you for raising this question. Although scRNA-seq and scATAC-seq\nare separate biological modalities, our model processes both in their peak-count\nrepresentations. We have adapted the embedding module of ChrisCell to encode\nchromatin regions, and we fine-tune the model by having it reconstruct the original\nscATAC-seq data, thus enabling it to learn from this modality.\n[1] Yu, Jixiang, et al. \"Unsupervised gene-cell collective representation learning with optimal\ntransport.\" Proceedings of the AAAI Conference on Artificial Intelligence. Vol. 38. No. 1. 2024.\n`"}, {"round": 3, "reviewer_comments": {"reviewer_1": "I thank the authors for their effort to revise the manuscript considering my comments. I have no further questions.", "reviewer_2": "Thank the authors for answering my questions.", "reviewer_3": "My comments have been well addressed.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": ""}]

Decision letter