Weighted Diversified Sampling for Efficient Data-Driven Single-Cell Gene-Gene Interaction Discovery

## Abstract

Abstract Gene-gene interactions play a crucial role in the manifestation of complex human diseases. Uncovering significant gene-gene interactions is a challenging task. Here, we present an innovative approach utilizing data-driven computational tools, leveraging an advanced Transformer model, to unearth noteworthy gene-gene interactions. Despite the efficacy of Transformer models, their parameter intensity presents a bottleneck in data ingestion, hindering data efficiency. To mitigate this, we introduce a novel weighted diversified sampling algorithm. This algorithm computes the diversity score of each data sample in just two passes of the dataset, facilitating efficient subset generation for interaction discovery. Our extensive experimentation demonstrates that by sampling a mere 1% of the single-cell dataset, we achieve performance comparable to that of utilizing the entire dataset.

## 1 Introduction

Gene-gene interactions play a crucial role in the manifestation of complex human diseases, including multiple sclerosis (Brassat et al., ; Motsinger et al., ; Slim et al., ), pre-eclampsia (Li et al., ; Diab et al., ; Williams and Pipkin, ; Oudejans and Van Dijk, ), and Alzheimer’s disease (Ghebranious et al., ; Hohman et al., ). Computational tools equipped with machine learning (ML) prove effective in uncovering these significant gene interactions (McKinney et al., ; Cui et al., ; Yuan and Bar-Joseph, ; Wei et al., ; Upstill-Goddard et al., ). By learning an ML model on massive single-cell transcriptomic data, we can identify gene-gene interactions associated with complex but common human diseases. Existing models rely on prior knowledge such as transcription factors (TF) (Wang et al., ; Yuan and Bar-Joseph, ; Chen et al., ; Shu et al., ) or existing gene-gene interaction (GGI) networks  (Ata et al., ; Yuan and Bar-Joseph, ), to infer new relationships. Although GGI networks and TFs are crucial for mapping biological processes, they frequently suffer from high false-positive rates and biases, particularly in large-scale in vitro experiments (Mahdavi and Lin, ; Rasmussen and et al., ). In response to these challenges, we propose that gene-gene interactions can be uncovered using purely data-driven methods.

The Rise of Transformers on Single-Cell Transcriptomic Data. Recent advances in natural language processing, particularly the development of Transformer models (Vaswani et al., ), have demonstrated significant potential in biological data analysis  (Hao et al., ; Theodoris et al., ; Bian et al., ; Cui et al., ). Transformer models are known for their ability to capture the dependencies between gene expressions. The information fused through the self-attention mechanism (Vaswani et al., ) is particularly suited for analyzing the intricate relationships in single-cell transcriptomic data. On the other hand, Transformer models also demonstrated superior performance when we scaled up their parameter size (Hao et al., ). This scaling capacity raises the researcher’s interest in training and deploying parameter-intensive Transformer models, denoted as single-cell foundation models (Cui et al., ). We would like to take this advantage for better gene-gene interaction discovery by identifying feature interactions within Transformer models.

Data-Driven Gene-Gene Interaction via Attention. In this work, we would like to advance the gene-gene interaction discovery with the Transformer models that have demonstrated superior performance on single-cell transcriptomic data. We see the self-attention mechanism  (Vaswani et al., ) as a pathway to facilitate the modeling of gene-gene interactions. In single-cell foundation models, the input to the model is a bag of $m$ gene expressions for a single cell. Next, in each layer and each head of the Transformer, there will be an attention map with shape $m\times m$ generated for this cell. Each entry of this attention map represents the interaction between two genes in this layer and this head. Assuming that we have a perfect Transformer that takes a cell gene expressions and correctly predicts if it is infected by a disease, we view the attention map of this cell as a strong indicator of disease-oriented gene-gene interactions.

Efficiency Challenge in Data Ingestion. Despite the transformative capabilities of Transformer models, one significant challenge remains: the efficient ingestion and processing of massive volumes of single-cell transcriptomic data. We are utilizing Transformer models with parameter sizes that exceed the hardware capacity, particularly that of the graphics processing unit (GPU). As a result, given a pre-trained Transformer, we have to perform batch-size computation on a massive single-cell transcriptomic dataset for computing gene-gene interactions through attention maps. This batch-size computation significantly enlarges the total execution time for scientific discovery. Moreover, the hardware in the real-world deployment environment for gene-gene interaction detection may have even more limited resources. Therefore, the current computational framework cannot support gene-gene interaction discovery on real-world single-cell transcriptomic datasets.

Our Proposal: Two-Pass Weighted Diversified Sampling.
In this paper, we introduce a novel weighted diversified sampling algorithm. This randomized algorithm computes the diversity score of each data sample in just two passes of the dataset. The proposed algorithm is highly memory-efficient and requires constant memory that is independent of the cell dataset size. Our theoretical analysis suggests that this diversity score estimates the density of the $\mathsf{Min}\text{-}\mathsf{Max}$ kernel defined on the cell-level gene expressions, which provides the foundation and justification of the proposed strategy. Through extensive experiments, we demonstrate how the proposed sampling algorithm facilitates efficient subset generation for interaction discovery. The results show that by sampling a mere 1% of the single-cell dataset, we can achieve performance comparable to that of utilizing the entire dataset.

Our Contributions. We summarize our contributions as fellows.

- •
We introduce a computational framework that advances the discovery of significant gene-gene interactions with CelluFormer, our proposed Transformer model that is trained on single-cell transcriptomic data.
- •
We pinpoint the challenge in data ingestion for the data-driven gene-gene interaction. Moreover, we argue that we should perform diversified sampling that selects a representative subset of single-cell transcriptomics data to fulfill the objective.
- •
We develop a diversity score for every cell in the dataset based on the $\mathsf{Min}\text{-}\mathsf{Max}$ kernel density. Moreover, we perform a randomized algorithm that efficiently estimates the $\mathsf{Min}\text{-}\mathsf{Max}$ kernel density for each cell. Furthermore, we use the estimated density to generate an effective subset for gene-gene interaction.

## 2 Data-Driven Single-Cell Gene-Gene Interaction Discovery

In this section, we propose a computing framework to perform gene-gene interaction discovery on single-cell transcriptomic data. We start by introducing the format of single-cell transcriptomic data. Next, we propose the formulation of our CelluFormer model tailored to single-cell data. Next, we present our multi-cell-type training to build an effective transformer model on single-cell data. Finally, given a pre-trained transformer, we showcase how to perform gene-gene interaction discovery by analyzing the attention maps.

### 2.1 Single-Cell Transcriptomic Data

Single-cell transcriptomic is a technology that profiles gene expression at the individual cell level. The profiled results, namely single-cell transcriptomic data, provide a unique landscape of gene expressions. In contrast to traditional bulk RNA-seq analysis, single-cell transcriptomic data allows for cell-level sequencing, which captures the variability between individual cells (Ata et al., ). Leveraging this high-resolution data allows scientists to gain insights into developmental processes, disease mechanisms, and cellular responses to environmental changes.

The single-cell transcriptomic data can be formulated as a set of high-dimensional and sparse feature vectors. We denote a single-cell transcriptomic dataset at $X$, where each cell $x\in X$ is a sparse vector with dimensionality $V\in\mathbb{N}_{+}$. Here $V$ represents the total number of genes we can observe in $X$. Since cell $x\in\mathbb{R}^{V}$ is a sparse vector, we can represent $x$ as a set $\{(i_{1},v_{1}),(i_{2},v_{2}),\cdots,(i_{k},v_{k})\}$. In this set, every tuple $(i,v)$ represents the expression of gene $i\in[V]$ with expression level $v\in\mathbb{R}$. Besides we can also denote cell $x$ as $[x_{1},x_{2},\cdots,x_{V}]$, where most of the $x_{i}$s are zeros.

In this data formulation, single-cell transcriptomic data for each cell is represented as a set of gene expressions, with different cells expressing varying genes. Additionally, even when two cells express the same gene, their expression levels may differ. Our research objective is to identify gene-gene interactions within the vocabulary $V$ that drive complex biological processes and disease mechanisms.

### 2.2 CelluFormer: A Single-Cell Transformer

Here, we propose our Transformer architecture, CelluFormer, to learn gene-gene interactions within single-cell transcriptomic data. Based on the set formulation of single-cell transcriptomic data, we believe that the order of genes is arbitrary and biologically meaningless.
Similar to scGPT (Cui et al., ), and scFoundation (Hao et al., ), our method adopts a permutation-invariant design. We define our permutation-invariant condition as follows.

###### Condition 2.1 .

Let $X$ denote a single-cell transcriptomic dataset. Given a single-cell data of cell $x\in X$, denoted as a set $\{(i_{1},v_{1}),(i_{2},v_{2}),\cdots,(i_{k},v_{k})\}$, a function $f:X\to\mathbb{R}$ should satisfy that, for any permutation $\pi$, $f(x)=f(\pi(x))$.

Figure: Figure 1: Distribution of Sequence Lengths in L6_CT Cell Type Data.
Refer to caption: extracted/5905729/figures/length_distribution_L6_CT.png

We see Condition 2.1 as a fundamental difference between the proposed Transformer and the sequence Transformers (Vaswani et al., ) widely used in natural language processing. For sequence Transformers, we have to ingest sequential masks during the training to ensure that the current token does not interact with the future token. Additionally, during the inference, the sequence Transformer should perform a step-by-step generation for each token. As a result, the sequence Transformer does not satisfy Condition 2.1. Moreover, the difference between CelluFormer and a vision Transformer (Dosovitskiy et al., ) is that the vision Transformer has a fixed sequence length for every input data sample. However, the number of genes expressed in each cell can vary a lot. For example, according to Figure 1, the number of genes expressed in a single cell can be up to 12,000 or more. Thus, we utilize a padding mask for the classification downstream task. Additional details regarding the implementation of CelluFormer are provided in Appendix C.1.

###### Condition 2.1 .

### 2.3 Multi-Cell-Type Training of CelluFormer

**Table 1: Performance comparison of models on neuronal cell dataset.**
| Model | Training Dataset | F1 Score | Accuracy |
| --- | --- | --- | --- |
| MLP | Pax6 | 78.91 | 82.71 |
| L5_ET | 62.02 | 73.31 |  |
| L6_CT | 91.14 | 92.01 |  |
| L6_IT_Car3 | 95.34 | 95.51 |  |
| L6b | 86.01 | 88.76 |  |
| Chandelier | 81.66 | 84.56 |  |
| L5_6_NP | 89.33 | 90.42 |  |
|  | All Neuronal Cell Types | 97.23 | 97.25 |
| CelluFormer | All Neuronal Cell Types | 98.12 | 98.12 |

We observe that there is a significant performance difference between Transformer models if we feed them with different styles of single-cell transcriptomic data. It is known that cells can be categorized into different types based on their functionality. For instance, neuronal cells represent the cell types that fire electric signals called action potentials across a neural network (Levitan and Kaczmarek, ). Our study suggests that Transformers should be trained on single-cell transcriptomic data from various cell types to achieve better performance. We showcase an example in Table 1. We train a Transformer model to classify whether a cell is an Alzheimer’s disease-infected cell or not. According to our study, CelluFormer proposed in Section 2.2 trained on neuronal cells outperforms traditional multilayer perceptron (MLP) with downstream training on a single cell type. However, we do not see this gap when we perform training of CelluFormer on a single cell type. As a result, we see that the Transformers generally prefer massive exposure to the single-cell transcriptomic data.

### 2.4 Gene-Gene Interaction Discovery via Attention Maps

In this paper, we would like to accomplish the following objective.

###### Objective 2.2 (Gene-gene interaction discovery) .

Let $X$ denote a single-cell transcriptomic dataset. Let $\mathcal{V}$ denote the genes expressed in at least one $x\in X$. Let $f:X\to\mathbb{R}$ denote a permutation invariant (see Condition 2.1) CelluFormer. $f$ can successfully predict whether any $x\in X$ is infected by disease $D$. We would like to find a gene-gene pair $(v_{1},v_{2})$ that contributes the most to $f$’s performance in $X$. Here $v_{1},v_{2}\in\mathcal{V}$.

Figure: Figure 2: Gene-gene interaction modeling with attention maps.
Refer to caption: extracted/5905729/figures/inference.png

Figure: Figure 3: Accumulating multiple cells’ average attention maps.
Refer to caption: extracted/5905729/figures/scatter_add.png

We see the self-attention mechanism of Transformers on a cell’s set style gene expressions as a pathway to model gene-gene interactions. CelluFormer takes a cell $x$’s gene expressions and produces an attention map $A_{i,j}\in\mathbb{R}^{m\times m}$ at encoder block $i$ and attention head $j$. Here $m$ represents the number of genes expressed in cell $x$. Since Transformer architecture uses the Softmax function to produce $A_{i,j}$, we can view the $p$th row of $A_{i,j}$ as the interaction between gene $p$ and all other genes in $x$. As a result, an attention map is a natural indicator of gene-gene interactions. Moreover,
if we have a perfect Transformer that takes a cell $x$ gene expressions and correctly predicts if it is infected by a disease, we view the attention map of this cell as an indicator of disease-oriented gene-gene interactions. Following this path, we propose a gene-gene interaction modeling approach as illustrated in Figure 2. For each cell $x$, we represent it as a set and generate a bag of embeddings from the gene embedding table. Next, we use the expression levels of each gene as a scaling factor for each gene’s embedding. Next, we take the average attention maps of all layers and all heads to obtain a gene-gene interaction map in this cell.

In Objective 2.2, we would like to see not only the gene-gene interactions just for cell $x$ but also the statistical evidence of how two genes interact in the dataset $X$. As a result, we propose to accumulate multiple cells’ averaged attention maps as illustrated in Figure 3. For $X$, we initialize $Z_{0}\in 0^{V\times V}$ matrix as the overall attention map before aggregation and $M_{0}\in 0^{V\times V}$ as the overall frequency dictionary before aggregation. Next, for each cell $x$ in the dataset, we remove its diagonal value in its averaged attention map as it represents self-interaction. Next, we perform scatter addition operations that merge $x$’s averaged attention map back to $Z_{0}$. We let $Z_{ij}$ add the interaction value of gene $v_{i}$ and $v_{j}$ in the average attention map of cell $x$ obtained in the Transformer model. Simultaneously, to eliminate the dataset bias of expressed genes, we count the number of appearances for each gene pair in the dataset. Once again, we perform scatter addition to record the counts back to $M_{0}$. This is done by updating $M_{0}$ through scatter addition, where $M_{ij}=M_{ij}+1$ for every occurrence of the gene pair $(v_{i},v_{j})$ in the dataset. Finally, we rank the off-diagonal values in $Z$ where $Z_{ij}\leftarrow\frac{Z_{ij}}{M_{ij}}$ to retrieve the top gene-gene interaction.

###### Objective 2.2 (Gene-gene interaction discovery) .

## 3 Weighted Diversified Sampling

In this section, we start by showcasing the data-efficiency problem when we use the trained CelluFormer for gene-gene interaction discovery. Following this, we define a diversity score for each cell in the dataset and propose a two-pass randomized algorithm to efficiently compute it. Lastly, we propose a weighted diversified sampling strategy on massive single-cell data.

### 3.1 Data-Intensive Computation for Gene-Gene Interaction Discovery

As illustrated in Section 2.4, once we have a pre-trained CelluFormer that can successfully predict whether a cell is infected by a disease or not with its gene expressions, we can perform gene-gene interaction discovery by passing massive cells into this model and get the accumulated attention map as Figure 3. However, this process requires data-intensive computation. For every cell in the dataset, we first need to compute the average attention map as illustrated in Figure 2. Next, we perform aggregations as shown in Figure 3. It is known that CelluFormer uses plenty of trainable parameters to achieve good disease infection classification performance. As a result, the computation complexity for generating a cell’s averaged attention map is expensive. Moreover, since the attention map for cell $x$ is $m\times m$, where $m$ is the number of genes expressed in $x$. According to Figure 1, we see that $m$ can be 12,000 or more. These giant attention maps consume the limited high bandwidth memory (HBM) in the graphics processing unit. Therefore, we have to perform batch-wise computation on a massive cell dataset for computing gene-gene interactions.
Moreover, given the scale of the dataset, any sampling algorithm with a runtime that grows exponentially with the dataset size is impractical.

Figure: Algorithm 1 Two-Pass Algorithm for Estimating $\mathsf{Min}\text{-}\mathsf{Max}$ Density

### 3.2 Two-Pass Randomized Algorithm for Computing 𝖬𝗂𝗇 ⁢ - ⁢ 𝖬𝖺𝗑 𝖬𝗂𝗇 - 𝖬𝖺𝗑 \mathsf{Min}\text{-}\mathsf{Max} sansserif_Min - sansserif_Max Density

In this work, we would like to address this data-efficiency challenge by raising and asking the following research question: Can we find a representative and small subset from the large cell dataset and still perform successful gene-gene interaction discovery? Moreover, we would like the procedure for finding this small subset as efficient as possible.

We would like to answer this question by proposing a diversity score of a cell in the dataset. To begin with, we introduce the $\mathsf{Min}\text{-}\mathsf{Max}$ similarity between two cell’s gene expressions.

###### Definition 3.1 ( $\mathsf{Min}\text{-}\mathsf{Max}$ Similarity) .

Given two cell’s gene expressions, denoted as $x,y\in\mathbb{R}^{V}$ (see Section 2.1), we define their $\mathsf{Min}\text{-}\mathsf{Max}$ similarity as:
$\mathsf{Min}\text{-}\mathsf{Max}(x,y)=\frac{\sum_{i}^{V}\min(x_{i},y_{i})}{
\sum_{i}^{V}\max(x_{i},y_{i})}$.

According to the definition, $\mathsf{Min}\text{-}\mathsf{Max}(x,y)\in[0,1]$. Higher $\mathsf{Min}\text{-}\mathsf{Max}$ means that two cell’s gene expressions are closer to each other. $\mathsf{Min}\text{-}\mathsf{Max}$ is widely viewed as a kernel (Li, ; Li et al., ; Li and Li, ) in statistical machine learning. In this paper, we would like to define a kernel density on top of the $\mathsf{Min}\text{-}\mathsf{Max}$ similarity.

###### Definition 3.2 ( $\mathsf{Min}\text{-}\mathsf{Max}$ Density) .

Given a cell dataset $X$, for every $q\in X$, we define its $\mathsf{Min}\text{-}\mathsf{Max}$ density as:
$\mathcal{K}(q)=\sum_{x\in X}\varphi(q,x)$,
where $\varphi(q,x):\mathbb{R}\to\mathbb{R}$ is a monotonic increasing function along with $\mathsf{Min}\text{-}\mathsf{Max}(q,x)$ similarity defined in Definition 3.1.

We view $\mathsf{Min}\text{-}\mathsf{Max}(q)$ density as an indicator of how diverse $q$ is in $X$. Smaller $\mathsf{Min}\text{-}\mathsf{Max}(q)$ means that all other $x\in X$ may be less similar to $q$, making $q$ a unique cell. On the other hand, higher $\mathsf{Min}\text{-}\mathsf{Max}(q)$ means that $X$ has some cells that have similar gene expressions with $q$, making $q$ less unique. However, to compute $\mathsf{Min}\text{-}\mathsf{Max}(q)$ for every $q\in X$ following Definition 3.2, we have to compute all pairwise $\mathsf{Min}\text{-}\mathsf{Max}(x,y)$ for any $x,y\in X$, which results in an unaffordable $O(n^{2}\mathbb{NNZ}(X))$ time complexity, where $n$ is the size of $X$ and $\mathbb{NNZ}(X)$ is the maximum possible number of genes expressed in a cell $x\in X$. To reduce this $n^{2}$ computation, we propose a randomized algorithm that takes advantage of 0-bit consistent weighted sampling (CWS) (Li, ) hash functions.

###### Definition 3.3 (0-bit Consistent Weighted Sampling Hash Functions (Li, 2015a ; Li et al., 2021 ) ) .

Let $\mathcal{H}$ denote a randomized hash function family.
If we pick a $h\in\mathcal{H}$ at random, for any two cell expressions $x,y\in\mathbb{R}^{V}$, we have $\Pr[h(x)=h(y)]=\mathsf{Min}\text{-}\mathsf{Max}(x,y)+o(1)$.
Here every $h\in\mathcal{H}$ is a hash function that maps any $x\in X$ to an integer in $[0,B)$. We denote $B$ as the hash range.

Here the $o(1)$ is a minor additive term with complex form. For simplicity, we refer the readers to (Li et al., ), Theorem 4.4 for more details.

This work presents an efficient randomized algorithm that estimates $\mathsf{Min}\text{-}\mathsf{Max}$ density $\mathcal{K}(q)$ (see Definition 3.2) for every $q\in X$. As showcased in Algorithm 1, we initialize an array $A$ with all values as zeros. Next, we conduct a pass over $X$. In this pass, for every $x\in X$, we compute its hash values after $R$ independent hash functions. Next, we increment $A_{r,h_{r}(x)}$ with $1$. After this pass, we take another pass at the dataset, for every $x\in X$, we take an average over the $A_{r,h_{r}(x)}$ and build a density score $w_{x}$.
We would like to highlight that Algorithm 1 requires only two linear scans of the dataset. The time complexity for this algorithm is $O(n\mathbb{NNZ}(X))$, which is linear to the dataset. Moreover, we show that Algorithm 1 produces an estimator to $\mathsf{Min}\text{-}\mathsf{Max}$ density.

###### Theorem 3.4 ( $\mathsf{Min}\text{-}\mathsf{Max}$ Density Estimator, informal version of Theorem B.1 ) .

Given a cell dataset $X$, for every $q\in X$, we compute $w_{q}$ following Algorithm 1. Next, we have
$\mathbb{E}[w_{q}]=\sum_{x\in X}(\mathsf{Min}\text{-}\mathsf{Max}(x,q)+o(1))$,
where $\mathsf{Min}\text{-}\mathsf{Max}$ is the $\mathsf{Min}\text{-}\mathsf{Max}$ similarity defined in Definition 3.1.
As a result, $w_{q}$ is an estimator for $\mathsf{Min}\text{-}\mathsf{Max}$ density $\mathcal{K}(q)$ defined in Definition 3.2 with $\varphi(q,x)=\mathsf{Min}\text{-}\mathsf{Max}(x,q)+o(1)$.

We provide the proof of Theorem 3.4 in the supplementary materials.

###### Definition 3.1 ( 𝖬𝗂𝗇 ⁢ - ⁢ 𝖬𝖺𝗑 𝖬𝗂𝗇 - 𝖬𝖺𝗑 \mathsf{Min}\text{-}\mathsf{Max} sansserif_Min - sansserif_Max Similarity) .

###### Definition 3.2 ( 𝖬𝗂𝗇 ⁢ - ⁢ 𝖬𝖺𝗑 𝖬𝗂𝗇 - 𝖬𝖺𝗑 \mathsf{Min}\text{-}\mathsf{Max} sansserif_Min - sansserif_Max Density) .

###### Definition 3.3 (0-bit Consistent Weighted Sampling Hash Functions (Li, 2015a ; Li et al., 2021 ) ) .

###### Theorem 3.4 ( 𝖬𝗂𝗇 ⁢ - ⁢ 𝖬𝖺𝗑 𝖬𝗂𝗇 - 𝖬𝖺𝗑 \mathsf{Min}\text{-}\mathsf{Max} sansserif_Min - sansserif_Max Density Estimator, informal version of Theorem B.1 ) .

### 3.3 Weighted Diversified Sampling with Inverse 𝖬𝗂𝗇 ⁢ - ⁢ 𝖬𝖺𝗑 𝖬𝗂𝗇 - 𝖬𝖺𝗑 \mathsf{Min}\text{-}\mathsf{Max} sansserif_Min - sansserif_Max Density

We propose to use the inverse form of $\mathsf{Min}\text{-}\mathsf{Max}$ density in Definition 3.2 as a score for diversity. We define it as normalized inverse $\mathsf{Min}\text{-}\mathsf{Max}$ density as below.

###### Definition 3.5 (Inverse $\mathsf{Min}\text{-}\mathsf{Max}$ Density (IMD)) .

Given a cell dataset $X$, for every $q\in X$, we define its normalized inverse $\mathsf{Min}\text{-}\mathsf{Max}$ density as
$\mathcal{I}(q)=\mathsf{Softmax}(1/\mathcal{K}(q))$,
where $\mathcal{K}(q)$ is the $\mathsf{Min}\text{-}\mathsf{Max}$ diversity for $q$ in Definition 3.2, $\mathsf{Softmax}$ is the softmax function that takes over all cells in $X$.

We view the IMD $\mathcal{I}(q)\in[0,1]$ as a monotonic increasing function for the diversity of $q$. Higher $\mathcal{I}(q)$ means that all other $x\in X$ may be less similar to $q$, making $q$ a unique cell. Moreover, IMD can be directly used as a sample probability to generate a representative subset of $X$ for Objective 2.2. Given $X$, we perform sampling without replacement to generate a subset $X_{\mathsf{sub}}\subset X$, where $x\in X$ has the sampling probability $\mathcal{I}(x)$. The advantages of sampling with IMD (see Definition 3.5) can be summarized as follows.

- •
The IMD $\mathcal{I}(q)$ can be an effective indicator for how diverse $q$ is in dataset $X$.
- •
Computing IMD is an efficient one-shot preprocessing process with just two linear scans of $X$ with time complexity $O(n\mathbb{NNZ}(X))$, where $n$ and $\mathbb{NNZ}(X)$ is defined in Section 3.2.
- •
The memory complexity of computing IMD is $O(RB)$, which can be viewed as constant since it is independent of $n$ and $\mathbb{NNZ}(X)$.

In the following section, we would like to evaluate the empirical performance of IMD in selecting a representative subset out of a massive single-cell transcriptomic dataset while still maintaining effective performance in data-driven gene-gene interaction discovery powered by CelluFormer.

###### Definition 3.6 (Estimated Interaction Score with WDS) .

Let $Z_{x}(v_{i},v_{j})$ denote the interaction value of gene $v_{i}$ and $v_{j}$ in the average attention map of cell $x$ obtained in the CelluFormer. For dataset $X$, we perform a sampling where each cell $x\in X$ is sampled with probability $\mathcal{I}(x)$ (see Definition 3.5) and get a subset $X_{s}$. Next, we define the estimated interaction score between gene $v_{i}$ and $v_{j}$ learned from $X$ as:

$$ $\displaystyle\widetilde{Z}(v_{i},v_{j})=\frac{\sum_{x\in X_{s}}Z_{x}(v_{i},v_{ j})\cdot\mathcal{I}(x)}{\sum_{x\in X_{s}}\mathcal{I}(x)},$ $$

where $\widetilde{Z}(v_{i},v_{j})$ is an unbiased estimator for the expectation of $Z(v_{i},v_{j})$ in distribution with density $\mathcal{I}(x)$. Formally,

$$ $\displaystyle\mathbb{E}[\widetilde{Z}(v_{i},v_{j})]=\mathbb{E}_{x\sim\mathcal{ I}(x)}[Z_{x}(v_{i},v_{j})],$ $$

$$ $\displaystyle\mathbf{Var}[\widetilde{Z}(v_{i},v_{j})]=~{}\frac{\sum_{x\in X_{s }}\mathcal{I}(x)^{2}}{(\sum_{x\in X_{s}}\mathcal{I}(x))^{2}}\mathbf{Var}_{x \sim\mathcal{I}(x)}[Z_{x}(v_{i},v_{j})].$ $$

###### Definition 3.5 (Inverse 𝖬𝗂𝗇 ⁢ - ⁢ 𝖬𝖺𝗑 𝖬𝗂𝗇 - 𝖬𝖺𝗑 \mathsf{Min}\text{-}\mathsf{Max} sansserif_Min - sansserif_Max Density (IMD)) .

###### Definition 3.6 (Estimated Interaction Score with WDS) .

## 4 Experiment

In this section, we want to validate the effectiveness of our gene-gene interaction pipeline as well as the two-pass diversified sampling algorithm 1. There are a few research questions we want to answer:

- •
RQ1: How does the proposed Transformer-based computing framework introduced in Section 2 perform in gene-gene interaction discovery?
- •
RQ2: How does the $\mathsf{Min}\text{-}\mathsf{Max}$ density estimated by two-pass diversified sampling Algorithm 1 characterize the diversity of a cell in the whole dataset? Is this estimated $\mathsf{Min}\text{-}\mathsf{Max}$ density useful?
- •
RQ3: How does the estimated $\mathsf{Min}\text{-}\mathsf{Max}$ density perform in improving data-efficiency of gene-gene interaction discovery? How is the quality of the subset sampled according to the estimated $\mathsf{Min}\text{-}\mathsf{Max}$ density?

### 4.1 Settings

Dataset: For the training dataset, we employ the Seattle Alzheimer’s Disease Brain Cell Atlas (SEA-AD) (Gabitto et al., ), which includes single nucleus RNA sequencing data of 36,601 genes (as 36,601 features) from 84 senior brain donors exhibiting varying degrees of Alzheimer’s Disease (AD) neuropathological changes.
By providing extensive cellular and genetic data, SEA-AD enables in-depth exploration of the cellular heterogeneity and gene expression profiles associated with AD. To facilitate a comparative analysis between AD-affected and non-AD brains, we select cells from 42 donors classified within the high-AD category and 9 donors from the non-AD category, based on their neuropathological profiles. This selection criterion ensures a robust comparison, aiding in the identification of gene-gene interactions linked to AD progression (Gabitto et al., ).
The dataset is comprehensively annotated, covering 1,240,908 cells across 24 distinct cell types. We selected 18 neuronal cell types as our final training dataset since we believe neuronal cells are more likely to reveal explainable gene-gene interactions that are related to Alzheimer’s Disease compared to non-neuronal cells. To better detect expression relationships among genes, we apply the Seurat Transformation Function (Stuart et al., ) to eliminate the problem of sequence depth difference.

Model: For the SEA-AD dataset, we designed a CelluFormer model as explained in 2.2 to predict labels indicative of Alzheimer’s disease conditions. Further details can be found in the Appendix C.1.

Baselines:
Our proposed algorithm leverages the attention maps of the Transformer models. Accordingly, we compare our method with three statistical methods, Pearson Correlation, CS-CORE, and Spearman’s Correlation (Freedman et al., ; Su et al., ; De Smet and Marchal, ). While these methods are widely adopted by biologists for gene co-expression analysis, gene co-expression values alone do not provide information about the relationship between gene pairs and Alzheimer’s Disease. To identify gene-gene interactions relevant to Alzheimer’s Disease, we apply these methods to subsets containing disease and non-disease cells respectively, and calculate their gene co-expression values. The difference in co-expression values between disease and non-disease cells is then used as a final score to rank the gene pairs. We also present more experiments in Appendix D.1 that demonstrate how Transformers aggregate data with varying labels.

Our baseline includes NID (Tsang et al., ), a traditional feature interpretation technique that extracts learned interactions from trained MLPs. NID identifies interacting features by detecting strongly weighted connections to a standard hidden unit in MLPs after training. We evaluated our CelluFormer model against the MLP model, with performance results presented in Table 1.

Additionally, to comprehensively evaluate RQ1, we utilized two existing single-cell large foundation models to assess our algorithm. Specifically, we fine-tuned two foundation models, scFoundation (Hao et al., ) and scGPT (Cui et al., ), to classify whether a cell is AD or non-AD (performance results are provided in Table 4). We then applied our gene-gene interaction discovery pipeline using the attention maps of these foundation models.

In the sampling experiments, we compare WDS with uniform sampling since none of them requires preprocessing time exponential to the dataset size.

Evaluation Metric:
For a comprehensive evaluation encompassing the entire ranked list of gene-gene interactions, we utilized the Kolmogorov-Smirnov test, which was facilitated by the GSEApy package (Fang et al., ) in Python. We select normalized enrichment score (NES) (Subramanian et al., ) as our evaluation metric.
The ground truth dataset is sourced from BioGRID and DisGenet (Oughtred et al., ; Piñero et al., ).
For our experiments, we extract a subset of DisGenet that includes genes associated with Alzheimer’s Disease. We then filter out genes in BioGRID that are not present in this DisGenet subset. Finally, we obtain a filtered BioGRID dataset containing only genes relevant to Alzheimer’s Disease. We provide more explanations about our evaluation metrics in Appendix C.2.

### 4.2 The Effectiveness of Transformers in Gene-Gene Interaction Discovery (RQ1)

To evaluate the effectiveness of our Transformer-based framework for gene-gene interaction discovery, we performed feature selection across seven different cell types used as inference datasets. Additionally, we used a dataset encompassing all neuronal cell types to assess the overall performance of various models. As shown in Table 2, Transformer-based methods, including CelluFormer, scGPT and scFoundation, significantly outperformed other baselines.

This result indicates that our proposed Transformer-based framework is more effective and stable at extracting general and global gene-gene interaction information. In addition, the foundation models, scGPT and scFoundation, achieved comparable performances with other baselines across some of the datasets.
We attribute this outcome to two main factors. Overfitting to Pretrained Knowledge: A foundation model, particularly a large one, might have learned very general or specific knowledge during its pretraining phase. When fine-tuning for a specific task, the model might overfit the preexisting knowledge, leading to poor generalization of the new task data. Mismatch Between Pretraining and Fine-Tuning Data: If the data distribution for fine-tuning is significantly different from the data on which the foundation model was trained, the model might struggle to adapt, resulting in worse performance. A model trained from scratch on the specific task data may perform better as it directly optimizes for that data distribution.

**Table 2: Performance comparison of models on neuronal cell data. To evaluate different models on datasets with varying sizes, we further select 7 neuronal cell types from all neuronal cell types. CelluFormer, scGPT, scFoundation, MLP, Pearson Correlation, Spearman’s Correlation, and CS-CORE were tested on 8 different datasets to obtain their gene pair rankings.**
| Dataset | CelluFormer | scFoundation | scGPT | NID | Pearson | CS-CORE | Spearman |
| --- | --- | --- | --- | --- | --- | --- | --- |
| L5_ET | 1.15 | 1.04 | 1.23 | 0.90 | 0.50 | 1.11 | 0.91 |
| L6_CT | 1.18 | 1.03 | 1.17 | 1.54 | -0.21 | 1.06 | 0.72 |
| Pax6 | 1.25 | 0.82 | 1.01 | 1.04 | 0.93 | 0.95 | 1.15 |
| L5_6_NP | 1.21 | 1.06 | 1.50 | 1.49 | 0.87 | 0.92 | 0.95 |
| L6b | 1.13 | 0.99 | 1.23 | 0.62 | 0.75 | 0.62 | 1.08 |
| Chandelier | 1.17 | 1.16 | 1.09 | 1.07 | 0.94 | 1.06 | 0.96 |
| L6_IT_Car3 | 1.22 | 0.90 | 0.66 | 1.19 | 0.59 | 1.08 | 0.86 |
| All neuron data | 1.17 | 1.02 | 0.99 | 0.86 | 1.01 | 1.06 | 1.04 |

### 4.3 Ablation Studies (RQ2 & RQ3)

We addressed these questions by comparing our weighted diversified sampling (WDS) method with uniform sampling across various sample sizes, ranging from 1% to 10% of the original dataset. We generated data subsets for each cell type using WDS and uniform sampling. We then applied our Transformer-based framework for feature selection at each sample size. Since CelluFormer consistently outperformed other baselines, we selected it as our base model. We repeated Each experiment five times and recorded the NES scores as the results. To evaluate the sampling methods, we calculated the average NES score across the five experiments. We also computed the Mean Square Error (MSE) between the NES scores from the sampling experiments and the ground truth derived from the entire dataset, as shown in Table 2. The evaluation results are presented in Table 3. We note that WDS consistently produced higher NES scores compared to uniform sampling. As the sample size increased, the NES scores from uniform sampling began to converge with the ground truth. In contrast, the NES scores from WDS consistently remained close to the ground truth, even at smaller sample sizes. The result indicates that while WDS offers a significant advantage in small samples by enabling the Transformer to capture a broader range of genetic interactions, its benefits diminish as more data becomes available. Moreover, we find that for some cell types, smaller samples of data outperformed larger samples of data on NES. This suggests that: (1) single-cell transcriptomic data may contain noises that affect gene-gene interaction discovery, and, (2) some complex gene-gene interaction patterns in single-cell transcriptomic data cannot be interpreted directly through attention maps. We also provide a detailed study on the choice of parameter $R$ in Algorithm 1 in Appendix D.2.

**Table 3: Evaluation Results for the transformer over sample data. For each cell type, we performed 8 groups of down-sampling regarding 4 different sample sizes and 2 sampling methods. We let the transformer conduct inferences over the sample data and generate results.**
|  |  | Mean of NES | MSE of NES |  |  |
| --- | --- | --- | --- | --- | --- |
| Dataset | Sample Size | Uniform | WDS | Uniform | WDS |
| L5_ET | 1% | 0.90 | 0.95 | 0.0127 | 0.0082 |
| 2% | 0.89 | 1.17 | 0.0131 | 0.0001 |  |
| 5% | 1.02 | 1.19 | 0.0036 | 0.0003 |  |
| 10% | 0.87 | 1.07 | 0.0158 | 0.0012 |  |
| L6_CT | 1% | 0.85 | 1.19 | 0.0207 | 0.0000 |
| 2% | 1.05 | 1.18 | 0.0030 | 4.30e-05 |  |
| 5% | 0.93 | 1.23 | 0.0122 | 0.0006 |  |
| 10% | 0.91 | 1.21 | 0.0136 | 0.0002 |  |
| Pax6 | 1% | 0.94 | 1.08 | 0.0184 | 0.0053 |
| 2% | 1.03 | 1.18 | 0.0098 | 0.0009 |  |
| 5% | 0.98 | 1.20 | 0.0139 | 0.0004 |  |
| 10% | 1.06 | 1.17 | 0.0072 | 0.0012 |  |
| L5_6_NP | 1% | 0.90 | 1.13 | 0.0192 | 0.0016 |
| 2% | 1.15 | 1.11 | 0.0009 | 0.0021 |  |
| 5% | 1.02 | 1.20 | 0.0076 | 4.54e-06 |  |
| 10% | 1.01 | 1.17 | 0.0080 | 0.0004 |  |
| L6b | 1% | 0.79 | 1.17 | 0.0226 | 0.0004 |
| 2% | 0.76 | 1.14 | 0.0266 | 0.0000 |  |
| 5% | 0.88 | 1.20 | 0.0121 | 0.0009 |  |
| 10% | 1.20 | 1.21 | 0.0010 | 0.0014 |  |
| L6_IT_Car3 | 1% | 0.78 | 1.20 | 0.0384 | 0.0001 |
| 2% | 0.87 | 1.15 | 0.0242 | 0.0011 |  |
| 5% | 0.97 | 1.17 | 0.0123 | 0.0006 |  |
| 10% | 0.97 | 1.18 | 0.0123 | 0.0003 |  |

## 5 Related Work

Single-Cell Transformer Models. Single-cell RNA sequencing (scRNA-seq), or single-cell transcriptomics, enables high-throughput insights into cellular systems, amassing extensive databases of transcriptional profiles across various cell types for the construction of foundational cellular models  (Hao et al., ). Recently, there has emerged a large number of transformer models pre-trained for single-cell RNA sequencing tasks, including scFoundation  (Hao et al., ), Geneformer  (Theodoris et al., ), scMulan  (Bian et al., ), scGPT  (Cui et al., ). These foundation models have gained a progressive understanding of gene expressions and can build meaningful gene encodings over limited transcriptomic data. Yet, the previous work did not pay attention to pairwise gene-gene interactions. In our work, we would like to highlight a fundamental functionality of single-cell foundation models: we must use these models to perform data-driven scientific discovery.

Randomized Algorithms for Efficient Kernel Density Estimation.
Kernel density estimation (KDE) is a fundamental task in both machine learning and statistics. It finds extensive use in real-world applications such as outlier detection (Luo and Shrivastava, ; Coleman et al., ) and genetic abundance analysis (Coleman et al., ). Recently, there has been a growing interest in applying hash-based estimators (HBE)(Charikar and Siminelakis, ; Backurs et al., ; Siminelakis et al., ; Coleman et al., ; Spring and Shrivastava, ) for KDE. HBEs leverage Locality Sensitive Hashing (LSH)(Indyk and Motwani, ; Datar et al., ; Li et al., ) functions, where the collision probability of two vectors under an LSH function is monotonic relative to their distance measure. This property allows HBE to perform efficient importance sampling using LSH functions and hash table-type data structures. Furthermore, (Liu et al., ) extend KDE algorithms as a sketch for the distribution. However, previous works have not considered LSH for weighted similarity. In this work, we focus on designing a new HBE that incorporates the $\mathsf{Min}\text{-}\mathsf{Max}$ similarity (Li, ), a weighted similarity measure.

## 6 Conclusion

Gene-gene interactions are pivotal in the development of complex human diseases, yet identifying these interactions remains a formidable challenge. In response, we have developed a pioneering approach that utilizes an advanced Transformer model to effectively reveal significant gene-gene interactions. Although Transformer models are highly effective, their extensive parameter requirements often impede efficient data processing. To overcome this limitation, we have introduced a weighted diversified sampling algorithm. This innovative algorithm efficiently calculates the diversity score of each data sample across just two passes of the dataset. With this method, we enable the rapid generation of optimized data subsets for interaction analysis. Our comprehensive experiments illustrate that by leveraging this method to sample a mere 1% of the single-cell dataset, we can achieve results that rival those obtained using the full dataset, significantly enhancing both efficiency and scalability.

## Appendix

## Appendix A More Related Work on Gene-Gene Interaction Discovery

In this section, we provide a more detailed review of the existing work on gene-gene interaction discovery. For gene-gene interaction network construction, genome-wide association studies (GWAS) are widely adopted by biologists to study gene associations using single-nucleotide polymorphisms (SNPs)  [Uffelmann et al., ]. However, GWAS have high computational costs and are simply based on direct genotype-phenotype associations instead of wired graph structure. To address this problem, many graphic models have emerged in recent years [Ata et al., ]. Network-based methods regard gene-gene interaction discovery as a task to construct a homogeneous graph among genes. For example, PRINCE [Vanunu et al., ] and VAVIEN [Erten et al., ] apply random work to predict new edges on existing protein-protein interaction (PPI) or gene-gene interaction (GGI) knowledge graphs. VGAE [Singh and Lio’, ] and GCAS [Rao et al., ] explore the potential to incorporate GNN and auto-encoder structure in the GGI network. In addition, existing work like DeepDRIM [Chen et al., ] and CNNC [Yuan and Bar-Joseph, ] successfully improve the construction of GGI through inferencing gene associations on scRNA sequencing data and known transcription factors (TF). While GGI networks and TFs are instrumental for mapping biological processes, they are often plagued by high false-positive rates and context-dependent inaccuracies, especially when derived from large-scale in vitro experiments [Mahdavi and Lin, , Rasmussen and et al., ]. Existing methods that rely on pre-established protein-protein interaction (PPI) or transcription factor (TF) networks are prone to bias because they tend to reinforce known interactions, making it difficult to objectively uncover novel gene-gene interactions. In contrast, our method circumvents this issue by directly discovering GGIs from scRNA-seq data without dependence on prior network knowledge.

## Appendix B Proofs of Theorem 3.4

###### Theorem B.1 ( $\mathsf{Min}\text{-}\mathsf{Max}$ Density Estimator, formal version of Theorem 3.4 ) .

Given a cell dataset $X$, for every $q\in X$, we compute $w_{q}$ following Algorithm 1. Next, we have

$$ $\displaystyle\mathbb{E}[w_{q}]=\sum_{x\in X}(\mathsf{Min}\text{-}\mathsf{Max}( x,q)+o(1)),$ $$

where $\mathsf{Min}\text{-}\mathsf{Max}$ is the $\mathsf{Min}\text{-}\mathsf{Max}$ similarity defined in Definition 3.1.
As a result, $w_{q}$ is an estimator for $\mathsf{Min}\text{-}\mathsf{Max}$ density $\mathcal{K}(q)$ defined in Definition 3.2 with $\varphi(q,x)=\mathsf{Min}\text{-}\mathsf{Max}(x,q)+o(1)$.

###### Proof.

According to Theorem 2 in [Coleman et al., ], the expectation of $w_{q}$ should be:

$$ $\displaystyle\mathbb{E}[w_{q}]=\sum_{x\in X}\Pr_{h\sim\mathcal{H}}[h(q)=h(x)]$ $$

According to Definition 3.3, we have

$$ $\displaystyle\Pr_{h\sim\mathcal{H}}[h(q)=h(x)]=\mathsf{Min}\text{-}\mathsf{Max }(x,q)+o(1)$ $$

.

As a result,

$$ $\displaystyle\mathbb{E}[w_{q}]=\sum_{x\in X}(\mathsf{Min}\text{-}\mathsf{Max}( x,q)+o(1))$ $$

Moreover, since $\mathsf{Min}\text{-}\mathsf{Max}(x,q)+o(1)$ is a monotonic increasing function of $\mathsf{Min}\text{-}\mathsf{Max}(x,q)$. We say that $w_{q}$ is an estimator for $\mathsf{Min}\text{-}\mathsf{Max}$ density $\mathcal{K}(q)$ defined in Definition 3.2 with $\varphi(q,x)=\mathsf{Min}\text{-}\mathsf{Max}(x,q)+o(1)$.
∎

###### Theorem B.1 ( 𝖬𝗂𝗇 ⁢ - ⁢ 𝖬𝖺𝗑 𝖬𝗂𝗇 - 𝖬𝖺𝗑 \mathsf{Min}\text{-}\mathsf{Max} sansserif_Min - sansserif_Max Density Estimator, formal version of Theorem 3.4 ) .

###### Proof.

## Appendix C Experiment Details

### C.1 Model Implementations

Transformer Configurations: In this work, we used the standard multi-head self-attention introduced in [Vaswani et al., ]. We do not see the potential of the proposed blocks in  [Lee et al., ] in our setting. Moreover, we perform padding on each batch of training and inference of single-cell data. Accordingly, we introduce a padding mask in the attention mechanism to avoid computation on the padded position. For each input sequence, we represent them as embedding by a lookup table that maps a vocabulary of 36,601 genes to 128-dimensional vectors. Subsequently, the embedded data passes through 4 transformer encoder blocks. Each encoder block features 8 attention heads, to capture complex, non-linear relationships within the data. Finally, the output is fed into a linear layer that classifies the data labels. Here the label for the cell can be disease-oriented, such as whether this cell is from an Alzheimer’s disease patient. We represent each input sequence by employing a lookup table that transforms a comprehensive vocabulary of 36,601 genes into 128-dimensional embedding vectors. These vectors are subsequently processed through a series of 4 Transformer encoder blocks. Each encoder block is equipped with 8 attention heads, a 512-dimensional feedforward layer, and a dropout layer in a ratio of 0.1. The processed outputs are then directed to a linear classification layer, which is tasked with predicting labels indicative of Alzheimer’s disease conditions. We adopted the Adam Optimization Algorithm to minimize the loss function  Kingma and Ba []. The model is trained under a learning rate of 1e-5 and the batch size of our data-loader is set as 128. The testing results for the transformer after 3 epochs of training are given in Table 1.

MLP Configurations: The MLP consists of 2 hidden layers, with 128 and 256 hidden units respectively. Each hidden layer is followed by a dropout and a Softplus module. The MLP is trained under a learning rate of 1e-4 and the batch size of our data-loader is set as 128. We adopted the Adam Optimization Algorithm to minimize the loss function  Kingma and Ba []. The testing results for the MLP after 80 epochs of training are given in Table 1.

**Table 4: Complete Performance comparison of models on neuronal cell data.**
| Model | Training Dataset | F1 Score | Accuracy |
| --- | --- | --- | --- |
| MLP | Pax6 | 78.91 | 82.71 |
| L5_ET | 62.02 | 73.31 |  |
| L6_CT | 91.14 | 92.01 |  |
| L6_IT_Car3 | 95.34 | 95.51 |  |
| L6b | 86.01 | 88.76 |  |
| Chandelier | 81.66 | 84.56 |  |
| L5_6_NP | 89.33 | 90.42 |  |
|  | All Neuronal Cell Types | 97.23 | 97.25 |
| CelluFormer | All Neuronal Cell Types | 98.12 | 98.12 |
| scGPT | All Neuronal Cell Types | 93.85 | 94.32 |
| scFoundation | All Neuronal Cell Types | 97.38 | 97.39 |

Fine-tuning configurations for scFoundatoin and scGPT:
For fine-tuning scGPT, we use an LR of 1e-4 and a batch size of 64. We utilize a step scheduler down to 90% of the original learning rate every 10 steps. The training process converges after 6 epochs. For scFoundation, we use an LR of 1e-4 and a batch size of 32. We fine-tune scFoundation for 10 epochs. The performances of scFoundation and scGPT on classifying disease cells are shown in Table 4.

Implementation and Computation Resources: Our codebase and workflow are implemented in PyTorch Paszke et al. []. We trained and tested our workflow on a server with 8 Nvidia Tesla V100 GPU and a 44-core/88-thread processor (Intel(R) Xeon(R) CPU E5-2699A v4 @ 2.40GHz).

### C.2 Evaluation Metrics

The normalized enrichment score (NES) is the main metric used to analyze gene set enrichment outcomes Subramanian et al. []. This score quantifies the extent of over-representation of a ground truth dataset at the top of the ranked list of gene-gene interactions. That is, the higher the better. We can calculate NES by starting at the top of the ranked list and moving through it, adjusting a running tally by increasing the score for each gene-gene interaction in the ground truth dataset and decreasing it for others based on each gene-gene interaction’s rank. This process continues until we evaluate the entire ranked list to identify the peak score, which is the enrichment score. The BioGRID Dataset provides human protein/genetic interactions. Specifically, BioGRID contributes $204,831$ protein/genetic interactions that help verify the enrichment of genuine biological interactions in a ranked list of gene-gene interactions. DisGenet contains 429,036 gene-disease associations (GDAs), connecting 17,381 genes to 15,093 diseases, disorders, and abnormal human phenotypes Oughtred et al. [], Piñero et al. [].

## Appendix D More Experments

### D.1 Contrastive Ranking

Here, we also explore alternative strategies for aggregating attention maps. While Pearson Correlation, Spearman’s Correlation, and CS-CORE themselves cannot capture the information between gene pairs the the target disease, we believe Transformers learn the difference among data with varying labels. Hence, we do not need to calculate the difference between attention maps aggregated on data with varying labels. However, given that the Transformer is trained to classify disease cells, we hypothesize that it likely assigns significant attention to specific gene pairs within disease cells. To evaluate this, we applied our pipeline to three distinct datasets. The experimental results summarized in Table 5 show that our pipeline achieves improved NES when both disease and non-disease cells are used as inputs. These findings suggest that the Transformer benefits from data both positive and negative labels to provide a more comprehensive understanding of features.

**Table 5: This experiment involves three groups. In the first group, the Transformer only takes the disease cells for inference. We directly evaluate the ranked list given by aggregated attention map across disease cells. In the second group, we calculate the aggregated attention maps on the disease cells and the non-disease cells respectively. The final attention map is obtained by subtracting these two attention maps. The third group is to aggregate attention maps across the whole dataset.**
| Strategy | L5_ET | L6_CT | Pax6 | L5_6_NP | L6b | Chandelier | L6_IT_Car3 |
| --- | --- | --- | --- | --- | --- | --- | --- |
| AD cells | 1.09 | 1.09 | 0.98 | 0.78 | 1.13 | 0.90 | 0.89 |
| AD cells - Non-AD cells | 1.08 | 0.89 | 1.05 | 0.76 | 0.82 | 0.65 | 1.39 |
| All cells | 1.15 | 1.18 | 1.25 | 1.21 | 1.13 | 1.17 | 1.22 |

### D.2 Empirical Study on Parameter R 𝑅 R italic_R in Algorithm 1

**Table 6: The Mean value of NES results across 5 experiments on L5_ET, L6_CT, and Pax6 cell type datasets.**
| Dataset | Sample Size | Mean of NES |  |  |  |
| --- | --- | --- | --- | --- | --- |
| Uniform | WDS with R=100 | WDS with R=200 | WDS with R=500 |  |  |
| L5_ET | 1% | 0.90 | 1.02 | 0.95 | 0.93 |
| 2% | 0.89 | 1.17 | 1.17 | 0.97 |  |
| 5% | 1.02 | 0.97 | 1.19 | 1.11 |  |
| 10% | 0.87 | 1.01 | 1.07 | 1.07 |  |
| L6_CT | 1% | 0.85 | 1.19 | 1.19 | 1.11 |
| 2% | 1.05 | 1.21 | 1.18 | 1.09 |  |
| 5% | 0.93 | 1.13 | 1.23 | 1.21 |  |
| 10% | 0.91 | 1.23 | 1.21 | 1.20 |  |
| Pax6 | 1% | 0.94 | 1.13 | 1.08 | 1.17 |
| 2% | 1.03 | 1.22 | 1.18 | 1.19 |  |
| 5% | 0.98 | 1.21 | 1.20 | 1.19 |  |
| 10% | 1.06 | 1.19 | 1.17 | 1.22 |  |

During our experiments on WDS, we observed that the value of $R$ (see Algorithm 1) has a noticeable impact on NES performance. In Table 6 and Table 7, we evaluate three different $R$ values ranging from 100 to 500. The results demonstrate that increasing $R$ leads to a significant decline in NES. Although WDS with smaller $R$ values yields relatively higher NES, it tends to diverge from the NES calculated on the entire dataset.

**Table 7: The MSE of NES results across 5 experiments on L5_ET, L6_CT, and Pax6 cell type datasets. The MSE values are calculated according to the results in Table 2.**
| Dataset | Sample Size | MSE of NES |  |  |  |
| --- | --- | --- | --- | --- | --- |
| Uniform | WDS with R=100 | WDS with R=200 | WDS with R=500 |  |  |
| L5_ET | 1% | 0.0636 | 0.0408 | 0.0178 | 0.0477 |
| 2% | 0.0653 | 0.0005 | 0.0004 | 0.0339 |  |
| 5% | 0.0181 | 0.0014 | 0.0310 | 0.0018 |  |
| 10% | 0.0790 | 0.0062 | 0.0192 | 0.0064 |  |
| L6_CT | 1% | 0.1033 | 0.0002 | 0.0002 | 0.0046 |
| 2% | 0.0151 | 0.0001 | 0.0014 | 0.0070 |  |
| 5% | 0.0610 | 0.0028 | 0.0025 | 0.0013 |  |
| 10% | 0.0681 | 0.0011 | 0.0031 | 0.0007 |  |
| Pax6 | 1% | 0.0920 | 0.0264 | 0.0135 | 0.0057 |
| 2% | 0.0488 | 0.0047 | 0.0006 | 0.0027 |  |
| 5% | 0.0695 | 0.0022 | 0.0015 | 0.0027 |  |
| 10% | 0.0362 | 0.0058 | 0.0028 | 0.0008 |  |

## Peer Reviews

confidence: 4
contribution: 2
initial_rating: 5
presentation: 3
questions: Please see the weaknesses.
soundness: 3
strengths: - The proposed method significantly reduces computational requirements without sacrificing accuracy. This approach addresses a key challenge in handling large-scale single-cell transcriptomic data.
- By leveraging Transformer models (CelluFormer) for gene-gene interaction discovery, the paper effectively adapts state-of-the-art NLP techniques to bioinformatics.
- The extensive experimental validation across multiple datasets and comparison with various baselines (e.g., Pearson Correlation, Spearman’s Correlation) provides empirical support for the proposed method’s effectiveness and robustness in data efficiency.
summary: This paper introduces a novel computational framework designed to discover gene-gene interactions linked to complex diseases through single-cell transcriptomic data. Utilizing a Transformer model named CelluFormer, the authors address the challenge of computational efficiency by implementing a weighted diversified sampling algorithm. This algorithm allows the selection of a representative data subset by calculating a diversity score based on the Min-Max density kernel.
weaknesses: - 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.
- 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 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.
- 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.

---

confidence: 2
contribution: 1
initial_rating: 3
presentation: 1
questions: Refer to the weaknesses for the questions that need to be addressed. 
More general ones:
- What are the related methods for this problem?
- What is the model proposed here?
soundness: 1
strengths: N/A
summary: The paper propose a method that is related to transformer architecture to do gene-gene interaction discovery in single cell data. The key idea of the method is a density-based sampling method to reduce the data size.

Edit: after the rebuttal, I still keep the scores.
weaknesses: The paper has major problems that prevent me from understanding the method itself and its relations to related methods such as transformer and previous work on the same problem.
- Problem setting is not clear: It is stated that data set X given, so what is a disease D is in the dataset? What need to be learnt? How many dimensions are there in X, |V| or m? 
- Model is not clear: what is "f"? The proposed model is not described anywhere in the text. Where is the interaction map in the model? What does it mean by "f" can successfully predict a disease infection? What does it mean by a gene pair that contribute the most to "f"?
- I think the paper simply refers the model to transformer architecture but the data here vectors!?!
- Conditions 2.1 on permutation invariant is totally misleading. "f" is defined for vectors, should we permute vector elements? 
- While the data is said to be sparse, there are up to ... 12000 expressed genes in a cell, most of them have 2000+ expressed genes?
- The definition of Min-Max density is clumsy since it can refer to its root: kernel density estimation.

---

confidence: 4
contribution: 1
initial_rating: 1
presentation: 2
questions: 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.

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.

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.

On the face of it, the fact that you can achieve the same performance from 1% of the data as from 100% can have two interpretations: you made fanastic use of the 1% to achieve performance comparable to 100%, or you made terrible use of the 100% and only achieved performance as good as using 1%.  How do we know which is the case here?

How are the cell type labels in Table 1 derived? More detail needs to be given to make this experiment reproducible.
soundness: 2
strengths: The problem of detecting gene-gene interactions directly from scRNA-seq data is important.
summary: This paper aims to use transformer model trained from scRNA-seq data to identify gene-gene interactions.  The approach involves combining the attention matrices over all the layers of the transformer and then evaluating whether the resulting aggregated attention values give higher weights to known pairs of interacting genes.  To make the approach computationally feasible, the authors also develop a sketching procedure to select a representation subset of cells.
weaknesses: The main idea here -- that you can infer gene-gene interactions by looking at the attention map in the transformer -- is pretty obvious.

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

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

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, and the paper fails to cite any existing methods that tackle this problem (e.g., Hie et al., Cell Systems 2019; Yang et al., ACM-BCB 2020; Yi & Stanley, bioRxiv, 2023; Hao et al. Nat Biotech 2024).  These methods should be employed and compared against.

More generally, there are many existing methods for detecting gene-gene interactions from scRNA-seq data.  Prominent examples include SCENIC (Aibar Nature Methods 2017), GRNBoost2 (Moerman Bioinformatics 2019), PIDC (Chan Cell Systems 2017), SCODE (Matsumoto Bioinformatics 2017), SCRIBE (Moerman Nature Communications 2019).  This large literature is not cited, and none of the methods therein is compared against.

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.

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.

Meta-review: This paper introduces a transformer model for identifying gene-gene interactions linked to complex diseases through single-cell transcriptomic data. 

The reviewers found strengths in the idea of adapting state-of-the-art NLP models to the problem of gene-gene interactions.

Significant concerns were notes by some reviewers. The review raised concerns about the claim that CelluFormer has advantages over traditional transformer-based methods. The review also notes concerns about the targeting of the paper and biological and/or methodological significance. Concerns were also raised about the specifics of the type of gene-gene interactions. It's important to carefully scope the paper to ensure the claims rest firmly on the foundation of the evidence. Clarifying the types of interactions is important to ensure this happens. Reviews also pointed out key literature that should be noted in the paper.

The authors and reviewers engaged in an extensive dialogue which was beneficial to the review process and to future improvements to the manuscript. Given the concerns and after a thorough review of the paper and conversation in the review process, the paper seems to be below the bar of acceptance at this time.


It seems it would be good for the authors to follow MPdj's suggestion about the literature here. 

Reviewer 8Hw had a productive dialogue with the authors and raised some important points about targeting of the paper and biological and/or methodological significance. I tend to agree with the reviewer and it seems the point is settled. Thanks for the discussion on that point. 

The reviewer also made a point about specificity of gene-gene interactions. It's important to carefully scope the paper to ensure the claims rest firmly on the foundation of the evidence. Clarifying the types of interactions is important to ensure this happens.

 I tend to agree that further methodological description and empirical evidence would be needed to support that claim. I also tend to agree with the sentiment that while "GNN- and SVM-based methods may not be directly applicable to single-cell RNA data due to differences in prior biological assumptions or data structure", it may be going too far at this stage of the literature to dismiss them when there is abundant recent work in that direction.

Decision: Reject