We thank the reviewer for the comments. Please see the following clarification.

[W1 - Evaluation of sampling algorithms - Our proposed sampling algorithm is a diversified sampling approach with the same sampling complexity as uniformly random sampling.]

The proposed sketching procedure is not compared to any existing methods.

We demonstrate in Algorithm 1 that our diversified sampling approach requires only a linear pass over the dataset, significantly reducing preprocessing time compared to other algorithms with exponential preprocessing complexity in terms of the dataset size $n$. Furthermore, once preprocessing is complete, our sampling procedure operates in $O(1)$ time, matching the efficiency of uniform sampling.

As a result, in this paper, we focuses on how to achieve diversified sampling under the constaint that the sampling procedures operates in $O(1)$ time and the preprocessing procedures operates in $O(nd)$. We do not see other sampling algorithms based on kernel density follow this constaint. So we showcase a comparison with uniform sampling.

[W2 - Utilization of Transformer models - We highlight the usage of Transformer models for interaction modelling on set of tokens.]

The proposed transformer model is insufficiently described, and the paper doesn't say how it differs from existing models.

We present two key features of Celluformer. First, the model operates on a set of gene expressions without assuming a sequential order. Consequently, Celluformer avoids introducing positional biases, such as positional encodings or causal masking, in its architecture. Second, our objective is not to use the Transformer model for performance prediction. Instead, we leverage its attention maps to infer novel gene-gene interactions, utilizing the attention mechanisms as a tool for discovery. This approach represents an empirical advancement, applying state-of-the-art Transformer architectures for scientific exploration. This is described on line 224-237 and Figure 1 in the manuscript.

[W3 - Study on the mentioned literature - We highlight the foundamental difference in the settings with these approaches with a study.]

One of the major contributions here is a method for finding representative subsets of cells from scRNA-seq data (Section 3). Unfortunately, this problem is already fairly well studied ...

We would like to clarify our relationship with the mentioned literature.

First, we would like to clarify Section 3. Section 3 is about developing a novel sampling method for data efficiency when performing gene-gene interaction discovery. To our best knowledge, this sampling approach is novel and has not been proposed or published by others.

Secondly, we appreciate the reviewer’s recommendation of some existing methods in gene-gene interactions. However, there are some fundamental differences between these methods and our method:

1. The goal of SCENIC, GRNBoost2, SCODE and SCRIBE is to discover gene regulatory network, with explicit intension of using transcription factor-target gene concept framework to model the data for gene regulation discovery. On the other hand, our method takes a data driven approach to identify gene-gene relationships without framing such relationship into any biological concepts like gene regulatory network. Even though gene regulation is an important gene-gene relationship from transcription profile, there could be other subtle signals of gene-gene interaction beyond gene regulation. Therefore, the scope and conceptual framework of our work is different from those works.
2. SCENIC, GRNBoost2, SCODE and SCRIBE takes traditional statistical or machine learning approaches to discover gene regulatory networks. On the other hand, our method is based on modern deep learning framework, the Transformer architecture, with specific emphasis on attention between features. We believe this attention focus aligns very well for the goal of discovering gene-gene interactions.

For these reasons, we feel like these applications are in a different arena from ours in terms of identifying gene-gene relationships. Therefore, we didn’t include these works. In the future work, we can explicitly focus on gene regulatory network conceptual framework and make head-to-head comparison with these algorithms.

[W4 - Line 161 - Please see our clarification.]

In line 161, you say that you train a Transformer model to classify whether a cell is an "Alzheimer's disease-infected cell or not." First, Alzheimer's is not an infection. But more importantly, there is no way to label individual cells as being affected by Alzheimer's or not. I am guessing that this sentence should say that you are labeling cells based on whether they come from an individual with Alzheimer's disease. How exactly this is done should be clarified.

Thank you very much for point this out. We have updated our statement regarding the disease, especially AD. We train a Transformer model to classify whether a cell is from an Alzheimer’s patient or a healthy individual. We revise the statement in line 67-69, line 161, line 172, line 217, line 236, line 249 and line 381 accordingly.

[W5 - Evaluation on gene-gene interaction - Please see our clarification.]

The evaluation of the gene-gene interactions is problematic. The approach amounts to filtering a large set of known gene-gene interactions to include only those interactions that are implicated in Alzheimer's disease. This ignores the fact that many genes not involved in Alzheimer's disease continue to function and interaction with one another in the cell. It's not actually clear to me whether the evaluation considers pairs of genes not involved in Alzheimer's. It seems, from the description, like these pairs are treated as non-interacting.

We would like to make it clear that the gene-gene relationships are inferred from a Transformer model trained to differentiate cells obtained from Alzheimer’s patients versus healthy controls. Because the model learns to differentiate Alzheimer’s disease from healthy controls, only genes and gene-gene relationships that are significantly different between these two states are identified/stands out. Even though many other genes function and interact, but they behave the same in both Alzheimer’s cells and healthy cells, therefore are not picked up by the Transformer model since there is no contribution from these genes in terms of differentiating Alzheimer’s and normal cells. We have revised our manuscript to clarify this in line 67-69, line 161, line 172, line 217, line 236, line 249 and line 381.

[Q1- Transformer for gene-gene interactions - We would like to highlight the imporatnce of gene-gene interaction and the limitation of existing pre-trained Transformers.]

Why develop a new transformer model if the goal is just to identify gene-gene interactions? It seems like it would be easier to work with some existing model.

We wish to emphasize the distinction in research focus between AI and biomedical studies. Our goal is to leverage our developed method to advance the understanding of human diseases, addressing the fundamental challenge of uncovering disease mechanisms in biomedical research. Specifically, we trained our Transformer model on a disease-specific dataset to identify gene-gene relationships pertinent to Alzheimer’s disease.

While we acknowledge the potential for applying our attention aggregation pipeline to other Transformer models, such as scGPT and scFoundation, as demonstrated in Table 2, we did not observe significant advantages from these foundation models for our specific task. Nonetheless, our computational framework remains generalizable, and we welcome the incorporation of additional pre-trained models for scRNA data in future explorations.

[Q2- Celluformer, scGPT and scFoundation - Our gene-gene interaction is generalizable to different transformer models. And we see the advantage of Celluformer directly trained on disease data.]

How does the proposed Transformer model differ from scGPT, GeneFormer, and scFoundation? On line 124 we are told that they are similar to CelluFormer, but I don't see any indication of how they differ.

Our work focuses on developing a novel attention aggregation pipeline while analyzing Alzheimer’s disease datasets to gain insights into disease mechanisms. While fine-tuning pre-trained models like scGPT is feasible, these models are trained on data unrelated to Alzheimer’s disease, which may result in gene-gene relationships that lack relevance to this specific condition. We don't see advantages of pre-trained models in our task in Table 2.

In contrast, by training our own Transformer-based model, Celluformer, exclusively on Alzheimer’s disease data, we ensure that the identified gene-gene relationships are directly tied to the disease. This targeted training avoids contamination from unrelated datasets, providing insights that are both precise and disease-specific.

Pre-trained models are typically trained on vast and diverse datasets, producing aggregate gene-gene relationships that represent generalized summaries. However, because Celluformer is trained on Alzheimer’s disease versus healthy datasets, the relationships it identifies are uniquely reflective of the disease, free from influence by unrelated data sources.

[Q3- Clarification on scatter addition - Sure!]

How do the "scatter addition operations" described in line 231 work? This entire paragraph would benefit from a more formal treatment, since I found the textual description very hard to follow.

As shown in line 224-236 and Figure 2, running inference on a cell $x$ on the transformer will produce a sparse attention matrix $A$, where the non-zero entries represents the average of attention scores between two genes expressed in cell $x$ when it go through the transformer. We take the average of all sparse attention matrices obtained from each cell to represent gene-gene interactions we modeled from the model on the whole dataset $X$.

[Q4- Justification of Sampling - Sure!]

On the face of it, the fact that you can achieve the same performance from 1% of the data a...

We appreciate the feedback provided. Our performance evaluation is centered on comparisons with existing methods, demonstrating that our approach significantly outperforms others when utilizing 100% of the data (see Table 2). Building on this benchmark, we further explore ways to enhance data efficiency within this framework. In particular, we emphasize that identifying a representative subset offers an effective strategy for leveraging large datasets, which constitutes the primary contribution of this work.

[Q5- Details of cell types - Sure!]

How are the cell type labels in Table 1 derived? More detail needs to be given to make this experiment reproducible.

The labels of the cell types are provided by data generator. Specifically, the labels are at two levels. First, cells are labeled by the sample origin, either from Alzheimer’s Disease patients or healthy control. Second, the transcription profiles of sequenced cells are analyzed by data generators and shared. Cells are clustered and assigned 24 different cell types based on clustering results and the cell type specific marker genes. Our analysis focused on 18 types of neuron cells as they are most relevant to the Alzheimer’s disease – a neural degenerative disease. We have revised the manuscript to clarify these points on line 398.

Dear Reviewer 8Hw9

Thanks again for helping review our paper! Since we are approaching the end of the author-reviewer discussion period, would you please check our response and our revised paper and see if there is anything we could add to address your concerns?

We appreciate the time and effort you spend reviewing our paper.

[Clarification on the interactions go into BioGrid - Please see the following clarificaiton]

Regarding gene regulatory networks, it is true that there are multiple types of gene-gene interactions. For example, the gene products (proteins) can interact with each other, in which case it is protein-protein interaction, which also has a large literature. It seems like by being inclusive about all types of interactions that you aim to discover, you end up not having to compare to any previous methods, which feels sort of odd to me. Maybe it would help to clarify exactly what kinds of interactions go into BioGrid.

I am glad that this reviewer agrees with our points. Indeed, BioGRID is an open‐access database resource that houses manually curated protein and genetic interactions from multiple species including yeast, worm, fly, mouse, and human. However, even all types of interaction there may not cover all the possible gene-gene relationships. On the other hand, our approach is purely data-driven, that is, any gene pair that contributes significantly to differentiating AD from non-AD cells will be captured. Our contribution is to develop a true discovery tool for any gene-gene relationships by combining large amounts of data from scRNA-Seq data and the amazing ability of deep learning frameworks like Transformer that can pick up subtle signals from large amounts of data. Our evaluation did emphasize whether our model, through a purely data-driven approach, successfully uncovers gene-gene interactions supported by BioGRID.

Regarding comparisons to other methods, we are not bypassing this step. Rather, the limited availability of directly comparable methods and literature makes a fair comparison challenging. Therefore, we benchmark our approach against alternative methodologies addressing the same problem, as presented in Table 2.

We thank the reviewer for the comments. Please see the following clarifications.

[W1 - Clarification on data settings.]

In our study, $D$ represents the disease of interest, such as Alzheimer's disease (AD). The goal is to develop a model that classifies whether a cell originates from an AD-origin individual or not. Subsequently, we use the model's attention mechanisms to uncover gene-gene interactions.

The SEA-AD dataset comprises 1,240,908 cells, meaning the dataset $X$ includes 1,240,908 elements. Each cell's gene expression data is drawn from a total of 36,601 genes, denoted as $V$. For each cell, the number of expressed genes, represented as $m$ varies from 2000 to 5000. Consequently, each element in $X$ is a sparse feature vector with $V = 36,601$ dimensions and up to $m$ non-zero values.

We have included a detail introduction of our data in line 116-120. For further details about our dataset and experimental setup, we refer reviewers to Sections 2.1 and 4.1 of the manuscript.

[W2 - Clarification on model settings.]

We use $f$ to represent the Transformer model used to classify whether a cell $x \in X$ originates from an individual with Alzheimer's Disease (AD) or from a healthy person. The self-attention mechanism in $f$ captures pairwise interactions between expressed genes in cell $x$. We hypothesize that if the model $f$ accurately performs the proposed classification task, gene pairs with high interaction values across different attention blocks in $f$ could serve as indicators of novel gene-gene interactions associated with AD, as such gene pair contributes significantly in differentiating AD cells from non-AD cells.

We have modified the manuscript to address this concerns on line 161, line 172.

[W3 & W5 - Clarification on Tokenziation of Gene Expressions.]

We clarify the tokenization technique used for representing gene expressions in a cell $x$, considering two perspectives:

1. $x$ is represented as a high-dimensional, sparse, non-binary vector, where each dimension corresponds to a gene. Non-zero entries indicate the genes expressed in this cell, and their values represent the respective expression levels.
2. Alternatively, $x$ can be represented as a set of tuples $(\text{id}, \text{value})$, where $\text{id}$ denotes the gene expressed in cell $x$, and $\text{value}$ represents the expression level.

Yes, the model in our work is a transformer architecture.

In our approach, we adopt the second perspective. A cell is represented as a set of tuples $(i, v_i)$ where $i$ is the ith gene and $v_i$ is the expression value of that gene. In the model , we utilize an embedding table to map each gene to a vector representation. For each cell, the retrieved vectors from the embedding table are scaled by their corresponding expression levels, resulting in contextualized gene representations specific to the cell.

Human genome has over 25,000 genes. However, because each cell has very little mRNA expressed for each gene, the detection rate of these expressed genes is low, ranging from 10-20% (2000-5000 genes) in most cells, but there could be few cells with higher number of genes measured than the range above. We have modified the manuscript on line 149 to clarify this point.

[W4 - Clarification on permutation invariance.]

We represent cells as a set of tuples (gene id and their expression value) rather than as vectors, introducing the concept of permutation invariance based on our interpretation of a cell as an unordered collection of tuples. This approach acknowledges that the sequential order of these tuples may vary due to differences in sequencing techniques. By treating the tuples as a set rather than a sequence, we ensure that the model avoids introducing positional bias into the representation.

This tuple-based representation, along with its implications for permutation invariance, is further clarified in line 116.

[W6 - Clarification on Min-Max Kernel Density.]

We have simplified our Min-Max density definitions in line 291. The Min-Max kernel stands out from general kernels due to its ability to be efficiently linearized using probabilistic hashing algorithms. This feature serves as the algorithmic foundation of our approach. By applying a low-cost computational procedure, we hash each element of the massive dataset once, enabling the use of count statistics for kernel density estimation.

Dear Reviewer QUgF

Thanks again for helping review our paper! Since we are approaching the end of the author-reviewer discussion period, would you please check our response and our revised paper and see if there is anything we could add to address your concerns?

We appreciate the time and effort you spend reviewing our paper and the suggestions you have provided!

We thank the reviewers for the comments and suggestions to improve our paper. Please see the following clarification.

[W1 - Model-specific patterns of Transformers - Our experiments do observe the impact of variations in transformer architectures!]

The paper does not thoroughly discuss potential limitations or biases in using Transformer attention maps for gene-gene interaction discovery, such as how model-specific patterns may impact biological interpretability or generalizability.

Transformer-based models with self-attention mechanisms excel at capturing pairwise interactions between genes via their attention maps. However, we do see that variations in Transformer architectures influence how gene-gene interactions are modeled. For example, unlike Celluformer that does not introduce any specific structure, scFoundation [1] employs the kernel-based approximation Transformer variant, Performer [2], in its decoder. Notably, using scFoundation's attention maps for gene-gene interactions results in a performance decline in GSEA analysis. We will include a detailed analysis in the updated version of the paper.

[W2 - Adaptivity of sampling methods - The exepected effective sample size is at $O(\log^2(n))$ scale.]

While the diversity score and sampling algorithm are well-motivated, the paper lacks detailed explanations regarding parameter sensitivity (e.g., choice of sample size) and the scalability of the method to even larger datasets or different cell types beyond the focus on Alzheimer’s Disease data.

The weighted diversified sampling approach introduced in this work builds on theoretical foundations related to local density estimation [3]. According to this analysis, the sample size required to approximate the target kernel density with a multiplicative error of $\epsilon$ and failure probabiltiy $\delta$ bounded by $O(\log^2(n)\cdot \log(1/\delta)/\epsilon^2)$, where $n$ represents the total number of elements in the dataset. We have included a discussion in lines 345–350 analyzing the sample size required to estimate the interaction score using the proposed weighted diversified sampling approach.

[W3 - Comparison with non-Transformer-based approaches - We have included multiple baselines with our available computational resources.]

The paper could benefit from a more comprehensive comparison with existing gene interaction discovery techniques, especially non-Transformer-based methods that might offer complementary insights or efficiency advantages.

Table 2 presents the evaluation of our approach alongside non-transformer-based methods, including NID, Pearson, CS-CORE, and Spearman. Additional baselines, such as locCSN and SpQN, were not included due to their high computational cost, which exceeds the resources outlined in Appendix C.1.

[W4 - Analysis of scGPT performance - We agree that there may be multiple reasons.]

For the scGPT model, there are cases where it performs better than the proposed method on specific datasets. Therefore, simply attributing the foundation model's lower performance to overfitting to pretrained knowledge or a mismatch between pretraining and fine-tuning data seems insufficient to support the claim that it underperforms compared to the proposed method.

We appreciate the reviewer's insightful assessment. The data handling approaches of foundation models, such as using rank instead of absolute expression values in scGPT, combined with the vast datasets used for training, make it challenging to isolate all factors contributing to the observed low performance. Potential influences may include differences in gene vocabulary and local minima arising from model training dynamics. In addition, variations in transformer model architecture could also contribute to the performance bias and variation as suggested by W1.

We have revised the claims and discussion in lines 430–444 to address potential reasons why foundation models have not demonstrated significant advancements in capturing gene-gene interactions. Future research could focus on a comprehensive evaluation of these factors to better understand their impact on model performance, particularly in identifying gene-gene interactions. Emphasis should be placed on the role of gene vocabulary, the effects of training dynamics and the impact of model variation.

[1] Large-scale foundation model on single-cell transcriptomics [2] Rethinking Attention with Performers [3] Local Density Estimation in High Dimensions