# MuSe-GNN: Learning Unified Gene Representation From Multimodal Biological
Graph Data

Tianyu Liu
Yale University
<EMAIL_ADDRESS>
Yuge Wang
Yale University
<EMAIL_ADDRESS>
Rex Ying
Yale University
<EMAIL_ADDRESS>
Hongyu Zhao*
Yale University
<EMAIL_ADDRESS>
*: Corresponding author

###### Abstract

Discovering genes with similar functions across diverse biomedical contexts
poses a significant challenge in gene representation learning due to data
heterogeneity. In this study, we resolve this problem by introducing a novel
model called Multimodal Similarity Learning Graph Neural Network, which
combines Multimodal Machine Learning and Deep Graph Neural Networks to learn
gene representations from single-cell sequencing and spatial transcriptomic
data. Leveraging 82 training datasets from 10 tissues, three sequencing
techniques, and three species, we create informative graph structures for
model training and gene representations generation, while incorporating
regularization with weighted similarity learning and contrastive learning to
learn cross-data gene-gene relationships. This novel design ensures that we
can offer gene representations containing functional similarity across
different contexts in a joint space. Comprehensive benchmarking analysis shows
our model’s capacity to effectively capture gene function similarity across
multiple modalities, outperforming state-of-the-art methods in gene
representation learning by up to 97.5%. Moreover, we employ bioinformatics
tools in conjunction with gene representations to uncover pathway enrichment,
regulation causal networks, and functions of disease-associated or dosage-
sensitive genes. Therefore, our model efficiently produces unified gene
representations for the analysis of gene functions, tissue functions,
diseases, and species evolution.

###  1 Introduction

Progress in biological technology has broadened the variety of biological
data, facilitating the examination of intricate biological systems. A prime
example of such technology is single-cell sequencing, which allows for the
comprehensive characterization of genetic information within individual cells
[38, 78]. This technology provides access to the full range of a cell’s
transcriptomics, epigenomics and proteomics information, including gene
expression (scRNA-seq) [46, 114], chromatin accessibility (scATAC-seq) [20,
17], methylation [60, 68] and anti-bodies [85]. By sequencing cells from the
same tissue at various time points, we can gain insight into patterns of
cellular activity over time [30]. Moreover, spatial information of cellular
activity represents an equally vital, additional dimension. Such data are
defined as spatial data. All of these data are known as multi-omics data, and
integrating multi-omics data for combined analysis poses a significant
challenge.

However, the traditional idea of multi-omics data integration known as using
cells as anchors [56, 87, 39, 12, 21, 76] is only partially suitable because
of the following two challenges. 1. Different omics data pose their own
challenges. For example, the unit of observation of the spatial transcriptomic
data is different from other single-cell data, as a single spatial location
contains mixed information (see the left part of Figure 1 a) from various
cells [23] and it is not appropriate to generate spatial clusters based on
gene expression similarity [16]. Current research [86] also indicates that
chromosome accessibility feature is not a powerful predictor for gene
expression at the cell level. 2. The vast data volume in atlas-level studies
challenges high-performance computing, risking out-of-memory or time-out
errors [58, 101]. With nearly 37.2 trillion cells in the human body [77],
comprehensive analysis is computationally infeasible. 3. Batch effects may
adversely impact analysis results [54] by introducing noise. Consequently, an
efficient and powerful method focusing on multi-omics and multi-tissue data
(referred to as multimodal biological data) analysis is urgently needed to
address these challenges.

Acknowledging the difficulties arising from the cell-oriented viewpoint, prior
work shifted focus to the gene perspective. Using gene sets as a summary of
expression profiles, based on natural selection during the species evolution,
may provide a more robust anchor [8]. Protein-coding genes are also thought to
interact with drugs [19], which are more relevant to diseases and drug
discovery. Gene2vec [24] is a method inspired by Word2vec [18], which learns
gene representations by generating skip-gram pairs from the co-expression
network. Recently, the Gene-based data Integration and ANalysis Technique
(GIANT) [16] has been developed, based on Node2vec [34] and OhmNet [119], to
learn gene representations from both single-cell and spatial datasets.
However, as shown in Figure 1 (b), significant functional clustering for genes
from different datasets but the same tissue was not observed based on the gene
embeddings from these two models, because these methods do not infer the
similarity of associated genes from different multimodal data. Additionally,
they did not offer metrics to quantitatively evaluate the performance of gene
embeddings compared to baseline models.

Here we introduce a novel model called Multimodal Similarity Learning Graph
Neural Network (MuSe-GNN) 111Codes of MuSe-GNN:
https://github.com/HelloWorldLTY/MuSe-GNN for multimodal biological data
integration from a gene-centric perspective. The overall workflow of MuSe-GNN
is depicted in Figure 1 (a). Figure 1 (b) shows MuSe-GNN’s superior ability to
learn functional similarity among genes across datasets by suitable model
structure and novel loss function, comparing to GIANT and Gene2vec. MuSe-GNN
utilizes weight-sharing Graph Neural Networks (GNNs) to encode genes from
different modalities into a shared space regularized by the similarity
learning strategy and the contrastive learning strategy. At the single-graph
level, the design of graph neural networks ensures that MuSe-GNN can learn the
neighbor information in each co-expression network, thus preserving the gene
function similarity. At the cross-data level, the similarity learning strategy
ensures that MuSe-GNN can integrate genes with similar functions into a common
group, while the contrastive learning strategy helps distinguish genes with
different functions. Furthermore, MuSe-GNN utilizes more robust co-expression
networks for training and applies dimensionality reduction [89] to high-
dimensional data [73].

Figure 1: The workflow of MuSe-GNN, the visualization of gene embeddings, and
the problems of two exitsing methods, GIANT and Gene2vec. (a) The process of
learning gene embeddings by MuSe-GNN. Here we highlight the difference between
single-cell data and spatial data, and the major applications of gene
embeddings. Each spot from single-cell data represents one cell, while each
spot from spatial data represents a mixture of cells. (b) UMAPs [63] for the
gene embeddings of MuSe-GNN, Gene2vec and GIANT. They are colored by datasets.
Based on UMAPs we can conclude that both Gene2vec and GIANT failed to learn
the gene similarity based on the datasets from the same tissue, while MuSe-GNN
produces meaningful embeddings agnostic to datasets.

To the best of our knowledge, this is the first paper in gene representation
learning that combines the Multimodal Machine Learning (MMML) concept [9, 27]
with deep GNNs [49, 81] design. Both approaches are prevalent in state-of-the-
art (SOTA) machine learning research, inspiring us to apply them to the joint
analysis of large-scale multimodal biological datasets 222Download links:
Appendix M.. As application examples, we first used the gene embeddings
generated by MuSe-GNN to investigate crucial biological processes and causal
networks for gene regulation in the human body. We also applied our model to
analyze COVID and cancer datasets, aiming to unveil potential disease
resistance mechanisms or complications based on specific differentially co-
expressed genes. Lastly, we used gene embeddings from MuSe-GNN to improve the
prediction accuracy for gene functions. Here genes are co-expressed means
genes are connected with the same edge.

Given the lack of explicit metrics to evaluate gene embeddings, we proposed
six metrics inspired by the batch effect correction problem [54] in single-
cell data analysis. We evaluated our model using real-world biological
datasets from one technique, and the benchmarking results demonstrated MuSe-
GNN’s significant improvement from 20.1% to 97.5% in comprehensive assessment.
To summarize the advantages of our model, MuSe-GNN addresses the outlined
problem about cross-data gene similarity learning and offers four major
contributions: 1. Providing an effective representation learning approach for
multi-structure biological data. 2. Integrating genes from different omics and
tissue data into a joint space while preserving biological information. 3.
Identifying co-located genes with similar functions. 4. Inferring specialized
causal networks of genes and the relationships between genes and biological
pathways or between genes and diseases.

###  2 Related work

Co-expression Network Analysis. While direct attempts at joint analysis of
gene functions across modalities are limited, there is relevant research in
identifying correlation networks based on a single dataset. WGCNA [52] is a
representative method that employs hierarchical clustering to identify gene
modules with shared functions. However, as an early tool, WGCNA has restricted
functionality. Its inference of co-expression networks necessitates more
rigorous methods, and it cannot be directly applied to the analysis of
multimodal biological data.

Network Based Biological Representation Learning. Apart from directly
generating gene co-expression networks, learning quantitative gene
representations in a lower dimensional space may better describe gene-gene
relationships and facilitate downstream analyses. Gene2vec generates gene
embeddings based on the co-expression network from a given database. However,
it disregards the expression profile information and is based on the old Gene
Expression Omnibus (GEO) up until 2017 [24]. GIANT leverages Node2vec [34] and
OhmNet [119] to learn gene representations from both single-cell and spatial
datasets by constructing graphs and hypergraphs for training. However, this
approach still over-compresses multimodal biological datasets by removing
expression profiles. Moreover, their co-expression networks created by Pearson
correlation have a high false positive rate [82]. Additionally, some of the
datasets used by GIANT are of low quality (see Appendix A). There are also
methods sharing a similar objective of learning embeddings from other
datasets. GSS [71] aims to learn a set of general representations for all
genes from bulk-seq datasets [62] using Principal Component Analysis (PCA) and
clustering analysis. However, it is for bulk-seq data and cannot be directly
applied to single-cell datasets from different tissues. Gemini [103] focuses
on integrating different protein function networks, with graph nodes
representing proteins rather than genes.

Graph Transformer. In the deep learning domain, Transformer [95, 100] is one
of the most successful models, leveraging seq2seq structure design, multi-head
self-attention design, and position encoding design. Many researchers have
sought to incorporate the advantages of Transformer to graph structure
learning. TransformConv [81] introduces a multi-head attention mechanism to
the graph version of supervised classification tasks and achieved significant
improvements. MAGNA [98] considers higher-level hop relationships in graph
data to enhance node classification capabilities. Graphomer [110] demonstrates
the positive impact of the Transformer structure on various tasks using data
from the Open Graph Benchmark Large-Scale Challenge [44], which is further
extended by GraphomerGD [111]. Recently, GPS [75] proposes a general Graph
Transformer (GT) model by considering additional encoding types. Transformer
architecture also contributes solutions to several biological questions [113].
scBERT [108] generates gene and cell embeddings using pre-training to improve
the accuracy of cell type annotations. The substantial impact of these efforts
highlights the crucial contribution of the Transformer architecture to graph
data learning.

###  3 Methods

In the following sections of introducing MuSe-GNN, we will elaborate on our
distinct approaches for graph construction that utilize multimodal biological
data, followed by an explanation of our weight-sharing network architecture
and the elements of the final loss function.

####  3.1 Preliminaries

GNN. GNNs aim to learn the graph representation of nodes (features) for data
with a graph structure. Modern GNNs iteratively update the representation of a
node by aggregating the representations of its $k$-order neighbors ($k\geq 1$)
and combining them with the current representation. As described in [110],
considering a graph $G=\\{V,E\\}$ with nodes $V=\\{v_{1},v_{2},...,v_{n}\\}$,
the AGGREGATE-COMBINE update for node $v_{i}$ is defined as:

$a_{i}^{(l+1)}=\operatorname{AGGREGATE}^{(l+1)}\left(\left\\{h_{j}^{(l)}:j\in\mathcal{N}\left(v_{i}\right)\right\\}\right);h_{i}^{(l+1)}=\operatorname{COMBINE}^{(l+1)}\left(h_{i}^{(l)},a_{i}^{(l+1)}\right),$
(1)

where $\mathcal{N}\left(v_{i}\right)$ represents the neighbors of node $v_{i}$
in the given graph, and $h_{i}^{(l)}$ and $h_{i}^{(l+1)}$ represent the node
representation before and after updating, respectively.

Problem Definition. We address the gene embeddings generation task by handling
multimodal biological datasets, denoted as
$\mathcal{D}=\left(\left\\{V_{i},E_{i}\right\\}\right)_{i=1}^{T}$. Our goal is
to construct a model $\mathcal{M}(\cdot,\theta)$, designed to yield gene
embeddings set
$\mathcal{E}=\\{e_{1},...,e_{T}\\}=\mathcal{M}(\mathcal{D},\theta)$. In this
context, $\mathcal{D}$ represents the input, $\theta$ represents the
parameters, and $\mathcal{E}$ represents the output. In other words, we aim to
harmonize gene information from diverse modalities within a unified projection
space, thereby generating consolidated gene representations.

####  3.2 Graph construction

Before constructing gene graphs, our first contribution involves the selection
of highly variable genes (HVGs) for each dataset. These HVGs constitute a
group of genes with high variance that can represent the biological functions
of given expression profiles. Moreover, considering co-expression networks is
important for gene representation learning because it allows us to
characterize gene-gene relationships. As sequencing depth, or the total counts
of each cell, often serves as a confounding factor in the co-expression
networks inference [11], we employ two unique methodologies, scTransform [35]
and CS-CORE [88], to process scRNA-seq and scATAC-seq data, thus creating gene
expression profiles and co-expression networks unaffected by sequencing depth.
For spatial transcriptomic data, our focus is on genes displaying spatial
expression patterns. We use SPARK-X [116] to identify such genes and then
apply scTransform and CS-CORE. For detailed algorithmic information regarding
these methods, please see Appendix B. Additionally, we demonstrate the
immunity of CS-CORE to batch effects when estimating gene-gene correlations in
Appendix C. In all our generated graphs (equivalent to co-expression
datasets), nodes represent genes and edges represent co-expression relation of
genes.

####  3.3 Cross-Graph Transformer

To capitalize on the strengths of the Transformer model during our training
process, we integrate a graph neural network featuring a multi-head self-
attention design [81], called TransformerConv, to incorporate co-expression
information and generate gene embeddings. Details of TransformerConv can be
found in Appendix B.4. Incorporating multimodal information can estimate more
accurate gene embeddings, supported by Appendix D. The cross-graph transformer
can efficiently learn gene embeddings containing gene functions across
different graphs, advocated by the comparison of different network structure
choices in Appendix E.1.

Weight sharing. Given the variability among multimodal biological datasets, we
employ a weight-sharing mechanism to ensure that our model learns shared
information across different graphs, representing a novel approach for
learning cross-graph relationships. We also highlight the importance of
weight-sharing design in Appendix E.1.

Figure 2: The overall model architecture and the design of loss functions for
MuSe-GNN. The color of nodes in the green block represents common/different
genes across two datasets. The brown block represents the network architecture
of MuSe-GNN, and the blue block represents different loss function components
of MuSe-GNN. The color gradients of the left two matrices represent different
gene expression levels.

Datasets & Modalities Graph Transformer. Drawing inspiration from the hard
weight-sharing procedure [15, 115], we not only employ dataset-specific Graph
Transformer (GT) layers $L_{1},L_{2},...,L_{n}$ for each graph
($G_{1},G_{2},...,G_{n}$) from the same modality $m$, but also connect all
these dataset-specific layers to a set of shared GT layers, denoted as
$D_{m}$. This design showcases our novel approach to incorporating weight-
sharing into the GT framework. The forward process of MuSe-GNN, given dataset
$i$ with network parameter $\theta_{*}$, is defined as follows:

$X_{i}^{\prime}=D_{m}(L_{i}(G_{i};\theta_{L_{i}});\theta_{D_{m}}).$ (2)

Datasets Decoder. Here we propose a dataset-specific decoder structure based
on Multi-layer Perceptrons (MLP). This decoder model is crucial in
reconstructing the co-expression relationships among different genes,
showcasing our inventive use of MLP for this purpose. Given a graph $G_{i}$
and its corresponding gene embedding $e_{i}$, the decoding process of MuSe-
GNN, with network parameter $\theta_{dec,i}$, is defined as follows:

$E_{rec}=\textsc{DecoderMlp}(e_{i}e_{i}^{T};\theta_{dec,i}),$ (3)

where $E_{rec}$ represents the reconstructed co-expression network.

####  3.4 Graph Reconstruction Loss ($\mathcal{L}_{\text{BCE}}$)

Within a single graph, we implement a loss function inspired by the Graph
Auto-encoder (GAE) [48]. This function is designed to preserve two key
aspects: 1. the similarity among genes sharing common functions, and 2. the
distinctions among genes with differing functions. This innovative use of a
GAE-inspired loss function constitutes a significant contribution to the
methodological design. For a graph $G_{i}=\\{V_{i},E_{i}\\}$, the loss
function for edge reconstruction is defined as:

$\displaystyle
e_{i}=\textsc{EncoderGNN}(\\{V_{i},E_{i}\\};\theta_{enc})=D_{m}(L_{i}(G_{i};\theta_{L_{i}});\theta_{D_{m}}),$
(4) $\displaystyle
E_{rec}=\textsc{DecoderMlp}(e_{i}e_{i}^{T};\theta_{dec,i});E^{\prime}_{rec}=\textsc{flatten}(E_{rec});E^{\prime}_{i}=\textsc{flatten}(E_{i}),$
$\displaystyle\mathcal{L}_{\text{BCE}}=-\frac{1}{|E^{\prime}_{rec}|}\sum_{t=1}^{|E^{\prime}_{rec}|}\left[E^{\prime}_{i}[t]\cdot\log
E^{\prime}_{rec}[t]+\left(1-E^{\prime}_{i}[t]\right)\cdot\log\left(1-E^{\prime}_{rec}[t]\right)\right],$

where EncoderGNN and DecoderMlp represent the encoder and decoder parts of our
model, respectively. $\mathcal{L}_{\text{BCE}}$ denotes the computation of
binary cross entropy loss (BCELoss) for the given input data. $E_{rec}$ is the
reconstructed adjacency matrix and $|E^{\prime}_{rec}|$ is the length of its
flatten version. Further justification of the model design can be found in
Appendix F.

####  3.5 Weighted-similarity Learning ($\mathcal{L}_{\text{Sim}}$)

To integrate shared biological information across graphs from disparate
datasets, we fuse the reconstruction loss of the input graph structure with a
cosine similarity learning loss. In this process, we treat common HVGs between
each pair of datasets as anchors. Our aim is to maximize the cosine similarity
($\textsc{CosSim}(\textbf{a},\textbf{b})=\frac{\textbf{a}\cdot\textbf{b}}{||\textbf{a}||_{2}\cdot||\textbf{b}||_{2}}$)
for each common HVG across two arbitrary datasets (represented by the yellow
blocks in Figure 2), utilizing their gene embeddings. However, in practice,
different common HVGs may have different levels of functional similarity in
two datasets, which is hard to be directly quantified. Consequently, we employ
the shared community score as an indirect measurement, which is incorporated
as a weight for the cosine similarity of different common HVG pairs within the
final loss function. Considering two graphs $G_{i}=\\{V_{i},E_{i}\\}$ and
$G_{j}=\\{V_{j},E_{j}\\}$, and one of their shared genes $g$, we identify the
co-expressed genes of the given gene $g$ in both datasets, denoted as $N_{ig}$
and $N_{jg}$. Thus, the weight $\lambda_{ijg}$ for gene $g$ can be expressed
as follows:

$\lambda_{ijg}=\frac{|N_{ig}\cap N_{jg}|}{|N_{ig}\cup N_{jg}|}.$ (5)

We can iterate over all shared genes from 1 to $n$ between these two graphs,
ultimately yielding a vector as
$\lambda_{ij}=[\lambda_{ij1},..,\lambda_{ijn}]$. This vector encapsulates the
community similarity between the two graphs across all common HVGs. We can
then modify the cosine similarity of various gene pairs by first multiplying
this vector with cosine similarities and then summing the resultant values
across all genes. The negation of the outcome is our final weighted similarity
loss, denoted as $\mathcal{L}_{\text{Sim}}$. The ablation test for weighed
similarity loss can be found in Appendix E.1. More detailed explanations and
evidence supporting this design in multimodal conditions can be found in
Appendix G.

####  3.6 Self-supervised Graph Contrastive Learning
($\mathcal{L}_{\text{InfoNCE}}$)

Specifically, when integrating multimodal biological data, we employ the
contrastive learning strategy [72] to ensure that functionally similar genes
are clustered together as closely as possible, while functionally different
genes are separated apart from each other. We utilize Information Noise
Contrastive Estimation (InfoNCE) as a part of our loss function to maximize
the mutual information between the anchor genes and genes with the same
functions. This loss is applicable to different genes in two arbitrary graphs
during the training process. In general, if we represent the embeddings of $N$
genes as $\mathrm{Gene}_{N}=\\{e_{1},...,e_{N}\\}$, InfoNCE is designed to
minimize:

$\mathcal{L}_{\text{InfoNCE}}=-\mathbb{E}\left[\log\frac{\exp\left(e\cdot
k_{+}/\tau\right)}{\sum_{i=0}^{K}\exp\left(e\cdot k_{i}/\tau\right)}\right],$
(6)

where samples $\\{k_{0},k_{1},k_{2}...\\}$ compose a set of gene embeddings
known as keys of one dictionary and $e$ is a query gene embedding. $k_{+}$
represents a positive sample of $e$ and $k_{i}$ denotes a negative sample of
$e$. Equation 6 can be interpreted as a log loss of a $(K+1)$-way Softmax
classifier, which attempts to classify $e$ as $k_{+}$. $\tau$ is a temperature
parameter. $\tau$ is set to 0.07 referenced in MoCo [41].

####  3.7 Final Loss Function

In summary, the training objective of MuSe-GNN for graph $i$ comprises three
components:

$\displaystyle\min_{e_{i},e_{j}}$ $\displaystyle\text{
}\mathcal{L}_{\text{BCE}}(\textsc{DecoderMlp}_{i}(e_{i}e_{i}^{T}),E_{i})-\mathbb{E}\left[\textsc{CosSim}(e_{i}[\mathrm{Common}_{ij}],e_{j}[\mathrm{Common}_{ij}])\lambda_{ij}^{T}\right]$
(7)
$\displaystyle+\lambda_{c}\mathcal{L}_{\text{InfoNCE}}(e_{i}[\mathrm{Diff}_{i}]\oplus
e_{j}[\mathrm{Diff}_{j}],e_{i}[\mathrm{Diff}_{\mathcal{N}(i)}]\oplus
e_{j}[\mathrm{Diff}_{\mathcal{N}(j)}]),$

where $\mathrm{Common}_{ij}$ denotes the index set for common HVGs,
$\mathrm{Diff}_{i}$ represents the index set for different HVGs in graph $i$,
and $\mathrm{Diff}_{\mathcal{N}(i)}$ indicates the index set for neighbors of
$\mathrm{Diff}_{i}$ in graph $i$. To expedite the training process and
conserve memory usage, we sample graph $j$ for each graph $i$ during model
training. We also employ multi-thread programming to accelerate the index set
extraction process. $\lambda_{c}$ is the weight for InfoNCE loss. All of the
components in MuSe-GNN are supported by ablation experiments in Appendix E.1.
Details of hyper-parameter tuning can be found in Appendix E.2.

###  4 Experiments

Algorithm 1 Multimodal Similarity Learning Graph Neural Network (MuSe-GNN)

1:Input: Model EncoderGNN, DecoderMlp; Dataset
$\mathcal{D}=\left(\left\\{V_{i},E_{i}\right\\}\right)_{i=1}^{T}$; Number of
epochs $K$; Neighbor overlap list $\lambda$; Weight coefficient for
contrastive learning $\lambda_{c}$;

2:Helper Functions: A function to find the genes with the same name
FindCommonGenes; A function to find the genes with different name
FindDiffGenes; A function to calculate cosine similarity CosSim; A function to
sample data from a graph collection Sample; A optimizer function to update
weights Adam.

3:Output: Model EncoderGNN

4:INIT: initialize all parameters.

5:for $s$ in $K$ steps do

6: for $\text{id},\\{V,E\\}$ in enumerate$(\mathcal{D})$ do

7:
$e\leftarrow\textsc{EncoderGNN}(\\{V,E\\};\theta_{enc});E_{rec}\leftarrow\textsc{DecoderMlp}(ee^{T};\theta_{dec})$

8: $\mathcal{L}\leftarrow\mathcal{L}_{\text{BCE}}(E_{rec},E)$

9:
$\text{id}_{\text{new}},\\{V_{\text{new}},E_{\text{new}}\\}\leftarrow\textsc{Sample}(G\backslash\\{V,E\\})$

10:
$e_{\text{new}}\leftarrow\textsc{EncoderGNN}(\\{V_{\text{new}},E_{\text{new}}\\};\theta_{enc})$

11:
$e_{\text{Common}},e^{\prime}_{\text{Common}}\leftarrow\textsc{FindCommonGenes}(e,e_{\text{new}})$

12:
$\mathcal{L}\leftarrow\mathcal{L}-\textsc{CosSim}(e_{\text{Common}},e^{\prime}_{\text{Common}})\lambda_{\text{id},\text{id}_{\text{new}}}^{T}$

13:
$e_{\text{Diff}},e^{\prime}_{\text{Diff}},e_{\text{Diff}_{N}},e^{\prime}_{\text{Diff}_{N}}\leftarrow\textsc{FindDiffGenes}(e,e_{\text{new}})$

14:
$\mathcal{L}\leftarrow\mathcal{L}+\lambda_{c}\mathcal{L}_{\text{InfoNCE}}(e_{\text{Diff}}\oplus
e_{\text{Diff}_{N}},e^{\prime}_{\text{Diff}}\oplus
e^{\prime}_{\text{Diff}_{N}})$

15:
$\theta_{enc}\leftarrow\textsc{Adam}(\mathcal{L},\theta_{enc});\theta_{dec}\leftarrow\textsc{Adam}(\mathcal{L},\theta_{dec})$

16:Return EncoderGNN

Datasets & Embeddings generation. Information on the different datasets used
for different experiments is included at the beginning of each paragraph. The
training algorithm of our model is outlined in Algorithm 1. We stored the gene
embeddings in the AnnData structure provided by Scanpy [104]. To identify
groups of genes with similar functions, we applied the Leiden clustering
algorithm [91] to the obtained gene embeddings. Details can be found in
Appendix E.4.

Evaluation metrics. For the benchmarking process, we used six metrics: edge
AUC (AUC) [74], common genes ASW (ASW) [59], common genes graph connectivity
(GC) [59], common genes iLISI (iLISI) [59], common genes ratio (CGR), and
neighbors overlap (NO) to provide a comprehensive comparison. Detailed
descriptions of these metrics can be found in Appendix H. We computed the
metrics for the different methods and calculated the average rank (Avg Rank
$\in[1,9]$). Moreover, to evaluate the model performance improvement, we
applied min-max scaling to every metric across different models and computed
the average score (Avg Score $\in[0,1]$).

Baselines. We selected eight models as competitors for MuSe-GNN. The first
group of methods stems from previous work on learning embeddings for
biomedical data, including Principal Component Analysis (PCA) [74] used by
GSS, Gene2vec, GIANT (the SOTA model for gene representation learning),
Weight-sharing Multi-layer Auto-encoder (WSMAE) used by [109] and scBERT (the
SOTA model with pre-training for cell type annotation). The second group of
methods comprises common unsupervised learning baseline models, including GAE,
VGAE [48] and Masked Auto-encoder (MAE) (the SOTA model for self-supervised
learning) [42]. MuSe-GNN, GIANT, GAE, VGAE, WSMAE and MAE have training
parameter sizes between 52.5 and 349 M, and all models are tuned to their best
performance. Details are shown in Appendix E.2.

Biological applications. For the pathway analysis, we used Gene Ontology
Enrichment Analysis (GOEA) [1] to identify specific biological pathways
enriched in distinct gene clusters with common functions. Moreover, we used
Ingenuity Pathway Analysis (IPA) [2] to extract biological information from
the genes within various clusters, including causal networks [50] and sets of
diseases and biological functions. Biological pathways refer to processes
identified based on the co-occurrence of genes within a particular cluster.
The causal network depicts the relationships between regulatory genes and
their target genes. Disease and biological function sets facilitate the
discovery of key processes and complications associated with specific
diseases. Using gene embeddings also improves the performance of models for
gene function prediction. To visualize gene embeddings in a low-dimensional
space, we utilized Uniform Manifold Approximation and Projection (UMAP) [63].

Table 1: Avg Score across different tissues. Standard deviations are reported in Appendix E.3. Methods | Heart | Lung | Liver | Kidney | Thymus | Spleen | Pancreas | Cerebrum | Cerebellum | PBMC
---|---|---|---|---|---|---|---|---|---|---
PCA | 0.52 | 0.48 | 0.56 | 0.47 | 0.56 | 0.60 | 0.51 | 0.62 | 0.53 | 0.51
Gene2vec | 0.40 | 0.37 | 0.33 | 0.29 | 0.21 | 0.31 | 0.24 | 0.27 | 0.31 | 0.19
GIANT | 0.50 | 0.40 | 0.33 | 0.38 | 0.58 | 0.33 | 0.56 | 0.29 | 0.28 | 0.28
WSMAE | 0.50 | 0.47 | 0.54 | 0.46 | 0.57 | 0.53 | 0.52 | 0.55 | 0.59 | 0.50
GAE | 0.61 | 0.45 | 0.58 | 0.40 | 0.56 | 0.58 | 0.52 | 0.56 | 0.60 | 0.54
VGAE | 0.64 | 0.32 | 0.33 | 0.38 | 0.56 | 0.31 | 0.33 | 0.41 | 0.33 | 0.47
MAE | 0.36 | 0.47 | 0.50 | 0.45 | 0.41 | 0.52 | 0.39 | 0.50 | 0.49 | 0.50
scBERT | 0.41 | 0.49 | 0.55 | 0.62 | 0.17 | 0.58 | 0.46 | 0.60 | 0.61 | 0.58
MuSeGNN | 0.77 | 0.96 | 0.92 | 0.89 | 0.89 | 0.94 | 0.80 | 0.95 | 0.90 | 0.92

####  4.1 Benchmarking Analysis

We executed each method 10 times by using the same setting of seeds to show
the statistical significance based on datasets across different tissues. The
performance comparison of nine gene embedding methods is presented in Tables 1
and 14. Based on these two tables, MuSe-GNN outperformed its competitors in
terms of both average ranks and average scores across all the tissues. Based
on Table 1, for major tissues, such as heart and lung, MuSe-GNN’s performance
was 20.1% higher than the second-best method and 97.5% higher than the second-
best method in heart and lung tissue, respectively. According to Appendix E.3,
MuSe-GNN’s stability was also demonstrated through various metrics by
comparing standard deviations, including AUC, GC, and NO. In contrast, methods
such as Gene2vec, GAE, VGAE, MAE and scBERT exhibited significant instability
in their evaluation results for kidney or thymus. Consequently, we concluded
that MuSe-GNN is the best performing model for learning gene representation
based on datasets from different tissues, making it applicable to learn gene
embeddings from diverse multimodal biological data.

####  4.2 Analysis of Gene Embeddings from Multimodal Biological Data

In Figure 3, we displayed the integration results for multimodal biological
data from Humans. Figure 3 (a) and (b) demonstrated that MuSe-GNN could
successfully integrate genes from different modalities into a co-embedded
space, allowing us to identify functional groups using the Leiden algorithm
shown in Figure 3 (c). Furthermore, Figure 3 (d) revealed that most of the
clusters in (c) were shared across different modalities. We also identified
three significant functional groups in Figure 3 (a): the nervous system
(predominantly composed of the cerebrum and cerebellum [26]), the
cardiovascular system (mainly composed of heart, lung, and kidney [29]), and
the immunology system (primarily consisting of spleen, liver, and peripheral
blood mononuclear cells (PBMC) [70]). All systems are important in regulating
the life activities of the body. We also uncovered a pre-epigenetics group
(mainly consisting of scATAC-seq data without imprinting, modification, or
editing), emphasizing the biological gap existing in multi-omics and the
importance of post-transcriptional regulation.

Using GOEA, we could identify significant pathways enriched by different co-
embedded gene clusters. For instance, Figure 3 (d) displayed the top 5
pathways in an immunology system cluster. The rank was calculated based on the
negative logarithm of the false discovery rate. Since all top pathways were
related to immunological defense and response, it further supported the
accuracy of our embeddings in representing gene functions. For our analysis of
shared transcription factors and major pathways across all the tissues, please
refer to Appendix I. For our analysis of multi-species gene embeddings, please
refer to Appendix J.

Figure 3: Gene representation learning results for multimodal biological data.
(a) represents the UMAPs of gene embeddings colored by tissue type and
highlighted by biological system. (b) represents the UMAPs of gene embeddings
colored by omics type. (c) represents the UMAPs of gene embeddings colored by
common function groups. (d) is a Sankey plot [80] to show the overlap of
different modalities in the same clusters. (e) shows the top5 pathways related
to the genes in the special cluster discovered by GOEA. The bubble plots in
this paper were created based on ggplot2 [102].  Figure 4: Gene embeddings
from COVID samples and healthy samples. (a) represents the UMAPs of gene
embeddings colored by functional groups. (b) represents the UMAPs of gene
embeddings colored by datasets. (c) represents the gene embeddings colored by
the conditions, and the red circle reflects the differential co-expression
genes. (d) shows the top6 pathways related to the genes in the special cluster
discovered by GOEA. (e) represents the causal network existing in the special
cluster discovered by IPA. (f) represents the top diseases & biological
functions discovered by IPA.

####  4.3 Analysis of Gene Embeddings for Diseases

We generated gene embeddings for human pancreas cells from samples with and
without SARS-CoV-2 infection, as depicted in Figure 4 (a) and (b). We
identified specific genes from COVID samples that did not align with control
samples, which piqued our interest. These genes, highlighted by a red circle
in Figure 4 (c), could be interpreted as differentially functional genes in
diseased cells.

We conducted GOEA for the genes of interest and discovered a close
relationship among these gene enrichment results and the top 5 pathways
associated with immune activity. These results are displayed in Figure 4 (d).
For the genes within our target cluster, we utilized IPA to identify the
Entrez name of these genes, and 90.3% (122/135) genes in our cluster are
related to immunoglobulin. We could also infer the causal relationship
existing in the gene regulatory activity of the immune system. For example,
Figure 4 (e) showed a causal network inferred by IPA based on our genes
cluster. PARP16, as an enzyme, can regulate ERN1 and EIF2AK3, and certain
pathways are also related to this causal network. Moreover, we also showed the
relation between the set of genes and Disease & Bio functions in Figure 4 (f).
We identified top related Diseases & Bio functions ranked by negative
logarithm of p-value, and all of these diseases could be interpreted as
complications that may arise from new coronavirus infection [6, 64, 79, 69,
33, 7, 83, 84]. Our extra analyses for lung cancer data can be found in
Appendix K.

####  4.4 Analysis of Gene Embeddings for Gene Function Prediction.

Here we intend to predict the dosage-sensitivity of genes related to genetic
diagnosis (as dosage-sensitive or not) [90]. We used MuSe-GNN to generate gene
embeddings for different datasets based on an unsupervised learning framework
and utilized the gene embeddings as training dataset to predict the function
of genes based on k-NN classifier. k-NN classifier is a very naive model and
can reflect the contribution of gene embeddings in the prediction task.

In this task, we evaluated the performance of MuSe-GNN based on the dataset
used in Geneformer [90, 31], comparing it to the prediction results based on
raw data or Geneformer (total supervised learning). As shown in Table 2, the
prediction accuracy based on gene emebddings from MuSe-GNN is the highest one.
Moreover, the performance of gene embeddings from MuSe-GNN is better than
Geneformer, which is a totally supervised learning model. Such finding proves
the advantages of MuSe-GNN in the application of gene function prediction
task. Further application analysis can be found in Appendix L.

Table 2: Accuracy for dosage-sensitivity prediction | MuSe-GNN (unsup) | Geneformer (sup) | Raw
---|---|---|---
Accuracy | 0.77$\pm$0.01 | 0.74$\pm$0.06 | 0.75$\pm$0.01

###  5 Conclusion

In this paper, we introduce MuSe-GNN, a model based on Multimodal Machine
Learning and Deep Graph Neural Networks, for learning gene embeddings from
multi-context sequencing profiles. Through experiments on various multimodal
biological datasets, we demonstrate that MuSe-GNN outperforms current gene
embedding learning models across different metrics and can effectively learn
the functional similarity of genes across tissues and techniques. Moreover, we
performed various biological analyses using the learned gene embeddings,
leveraging the capabilities of GOEA and IPA, such as identifying significant
pathways, detecting diseases, and inferring causal networks. Our model can
also contribute to the study of the pathogenic mechanisms of diseases like
COVID and lung cancer, and improve the prediction performance for gene
functions. Overall, the gene representations learned by MuSe-GNN are highly
versatile and can be applied to different analysis frameworks.

At present, MuSe-GNN does not accept graphs with nodes other than genes as
input. In the future, we plan to explore more efficient approaches for
training large models related to Multimodal Machine Learning and extend MuSe-
GNN to a more general version capable of handling a broader range of
multimodal biological data.

###  6 Acknowledgements

This research was supported in part by NIH R01 GM134005, R56 AG074015, and NSF
grant DMS 1902903 to H.Z. We appreciate the comments, feedback, and
suggestions from Chang Su, Zichun Xu, Xinning Shan, Yuhan Xie, Mingze Dong,
and Maria Brbic.

### References

  * [1] Gene ontology enrichment analysis.  URL http://geneontology.org.
  * [2] Ingenuity pathway analysis.  URL https://www.qiagenbioinformatics.com/products/ingenuitypathway-analysis.
  * 10x Genomics [2018] 10x Genomics.  10x genomics acquires spatial transcriptomics.  _Science_ , 2018.  URL https://www.10xgenomics.com/news/10x-genomics-acquires-spatial-transcriptomics.
  * Altenhoff et al. [2021] Adrian M Altenhoff, Clément-Marie Train, Kimberly J Gilbert, Ishita Mediratta, Tarcisio Mendes de Farias, David Moi, Yannis Nevers, Hale-Seda Radoykova, Victor Rossier, Alex Warwick Vesztrocy, et al.  Oma orthology in 2021: website overhaul, conserved isoforms, ancestral gene order and more.  _Nucleic acids research_ , 49(D1):D373–D379, 2021.
  * Andrews and Hemberg [2018] Tallulah S Andrews and Martin Hemberg.  False signals induced by single-cell imputation.  _F1000Research_ , 7, 2018.
  * Aram et al. [2021] Khashayar Aram, Anant Patil, Mohamad Goldust, and Fateme Rajabi.  Covid-19 and exacerbation of dermatological diseases: a review of the available literature.  _Dermatologic Therapy_ , 34(6):e15113, 2021.
  * Ardestani Zadeh and Arab [2021] Arash Ardestani Zadeh and Davood Arab.  Covid-19 and male reproductive system: pathogenic features and possible mechanisms.  _Journal of molecular histology_ , 52:869–878, 2021.
  * Argelaguet et al. [2021] Ricard Argelaguet, Anna SE Cuomo, Oliver Stegle, and John C Marioni.  Computational principles and challenges in single-cell data integration.  _Nature biotechnology_ , 39(10):1202–1215, 2021.
  * Baltrušaitis et al. [2018] Tadas Baltrušaitis, Chaitanya Ahuja, and Louis-Philippe Morency.  Multimodal machine learning: A survey and taxonomy.  _IEEE transactions on pattern analysis and machine intelligence_ , 41(2):423–443, 2018.
  * Bartlett and Mendelson [2002] Peter L Bartlett and Shahar Mendelson.  Rademacher and gaussian complexities: Risk bounds and structural results.  _Journal of Machine Learning Research_ , 3(Nov):463–482, 2002.
  * Booeshaghi et al. [2022] A Sina Booeshaghi, Ingileif B Hallgrímsdóttir, Ángel Gálvez-Merchán, and Lior Pachter.  Depth normalization for single-cell genomics count data.  _bioRxiv_ , pages 2022–05, 2022.
  * Butler et al. [2018] Andrew Butler, Paul Hoffman, Peter Smibert, Efthymia Papalexi, and Rahul Satija.  Integrating single-cell transcriptomic data across different conditions, technologies, and species.  _Nature biotechnology_ , 36(5):411–420, 2018.
  * Cai et al. [2021] Tianle Cai, Shengjie Luo, Keyulu Xu, Di He, Tie-yan Liu, and Liwei Wang.  Graphnorm: A principled approach to accelerating graph neural network training.  In _International Conference on Machine Learning_ , pages 1204–1215. PMLR, 2021.
  * Cao et al. [2020] Junyue Cao, Diana R O’Day, Hannah A Pliner, Paul D Kingsley, Mei Deng, Riza M Daza, Michael A Zager, Kimberly A Aldinger, Ronnie Blecher-Gonen, Fan Zhang, et al.  A human cell atlas of fetal gene expression.  _Science_ , 370(6518):eaba7721, 2020.
  * Caruana [1993] R Caruana.  Multitask learning: A knowledge-based source of inductive bias1.  In _Proceedings of the Tenth International Conference on Machine Learning_ , pages 41–48. Citeseer, 1993.
  * Chen et al. [2022] Hao Chen, Nam D Nguyen, Matthew Ruffalo, and Ziv Bar-Joseph.  A unified analysis of atlas single cell data.  _bioRxiv_ , 2022.
  * Chen et al. [2018] Xi Chen, Ricardo J Miragaia, Kedar Nath Natarajan, and Sarah A Teichmann.  A rapid and robust method for single cell chromatin accessibility profiling.  _Nature communications_ , 9(1):1–9, 2018.
  * Church [2017] Kenneth Ward Church.  Word2vec.  _Natural Language Engineering_ , 23(1):155–162, 2017.
  * Cotto et al. [2018] Kelsy C Cotto, Alex H Wagner, Yang-Yang Feng, Susanna Kiwala, Adam C Coffman, Gregory Spies, Alex Wollam, Nicholas C Spies, Obi L Griffith, and Malachi Griffith.  Dgidb 3.0: a redesign and expansion of the drug–gene interaction database.  _Nucleic acids research_ , 46(D1):D1068–D1073, 2018.
  * Cusanovich et al. [2015] Darren A Cusanovich, Riza Daza, Andrew Adey, Hannah A Pliner, Lena Christiansen, Kevin L Gunderson, Frank J Steemers, Cole Trapnell, and Jay Shendure.  Multiplex single-cell profiling of chromatin accessibility by combinatorial cellular indexing.  _Science_ , 348(6237):910–914, 2015.
  * Domínguez Conde et al. [2022] C Domínguez Conde, C Xu, LB Jarvis, DB Rainbow, SB Wells, T Gomes, SK Howlett, O Suchanek, K Polanski, HW King, et al.  Cross-tissue immune cell analysis reveals tissue-specific features in humans.  _Science_ , 376(6594):eabl5197, 2022.
  * Dong and Kluger [2023] Mingze Dong and Yuval Kluger.  Towards understanding and reducing graph structural noise for gnns.  2023\.
  * Dong and Yuan [2021] Rui Dong and Guo-Cheng Yuan.  Spatialdwls: accurate deconvolution of spatial transcriptomic data.  _Genome biology_ , 22(1):145, 2021.
  * Du et al. [2019] Jingcheng Du, Peilin Jia, Yulin Dai, Cui Tao, Zhongming Zhao, and Degui Zhi.  Gene2vec: distributed representation of genes based on co-expression.  _BMC genomics_ , 20(1):7–15, 2019.
  * Du et al. [2021] Simon Shaolei Du, Wei Hu, Sham M. Kakade, Jason D. Lee, and Qi Lei.  Few-shot learning via learning the representation, provably.  2021\.  URL https://openreview.net/forum?id=pW2Q2xLwIMD.
  * Durmaz et al. [2001] Ramazan Durmaz, Metin Ant Atasoy, Gül Durmaz, Baki Adapinar, Ali Arslantaş, Aydin Aydinli, and Esref Tel.  Multiple nocardial abscesses of cerebrum, cerebellum and spinal cord, causing quadriplegia.  _Clinical neurology and neurosurgery_ , 103(1):59–62, 2001.
  * Ektefaie et al. [2023] Yasha Ektefaie, George Dasoulas, Ayush Noori, Maha Farhat, and Marinka Zitnik.  Multimodal learning with graphs.  _Nature Machine Intelligence_ , pages 1–11, 2023.
  * et al. [2019] HuBMAP Consortium et al.  The human body at cellular resolution: the nih human biomolecular atlas program.  _Nature_ , 574(7777):187–192, 2019.
  * Fessler [1997] HE Fessler.  Heart-lung interactions: applications in the critically ill.  _European Respiratory Journal_ , 10(1):226–237, 1997.
  * Fischer et al. [2019] David S Fischer, Anna K Fiedler, Eric M Kernfeld, Ryan MJ Genga, Aimée Bastidas-Ponce, Mostafa Bakhti, Heiko Lickert, Jan Hasenauer, Rene Maehr, and Fabian J Theis.  Inferring population dynamics from single-cell rna-sequencing time series data.  _Nature biotechnology_ , 37(4):461–468, 2019.
  * Franzén et al. [2019] Oscar Franzén, Li-Ming Gan, and Johan LM Björkegren.  Panglaodb: a web server for exploration of mouse and human single-cell rna sequencing data.  _Database_ , 2019:baz046, 2019.
  * Gabaldón and Koonin [2013] Toni Gabaldón and Eugene V Koonin.  Functional and evolutionary implications of gene orthology.  _Nature Reviews Genetics_ , 14(5):360–366, 2013.
  * Gayen Nee’Betal et al. [2022] Suhita Gayen Nee’Betal, Pedro Urday, Huda B Al-Kouatly, Kolawole Solarin, Joanna SY Chan, Sankar Addya, Rupsa C Boelig, and Zubair H Aghai.  Covid-19 infection during pregnancy induces differential gene expression in human cord blood cells from term neonates.  _Frontiers in Pediatrics_ , 10:547, 2022.
  * Grover and Leskovec [2016] Aditya Grover and Jure Leskovec.  node2vec: Scalable feature learning for networks.  In _Proceedings of the 22nd ACM SIGKDD international conference on Knowledge discovery and data mining_ , pages 855–864, 2016.
  * Hafemeister and Satija [2019] Christoph Hafemeister and Rahul Satija.  Normalization and variance stabilization of single-cell rna-seq data using regularized negative binomial regression.  _Genome biology_ , 20(1):1–15, 2019.
  * Hamilton et al. [2017] Will Hamilton, Zhitao Ying, and Jure Leskovec.  Inductive representation learning on large graphs.  _Advances in neural information processing systems_ , 30, 2017.
  * Han et al. [2018] Heonjong Han, Jae-Won Cho, Sangyoung Lee, Ayoung Yun, Hyojin Kim, Dasom Bae, Sunmo Yang, Chan Yeong Kim, Muyoung Lee, Eunbeen Kim, et al.  Trrust v2: an expanded reference database of human and mouse transcriptional regulatory interactions.  _Nucleic acids research_ , 46(D1):D380–D386, 2018.
  * Han et al. [2020] Xiaoping Han, Ziming Zhou, Lijiang Fei, Huiyu Sun, Renying Wang, Yao Chen, Haide Chen, Jingjing Wang, Huanna Tang, Wenhao Ge, et al.  Construction of a human cell landscape at single-cell level.  _Nature_ , 581(7808):303–309, 2020.
  * Hao et al. [2021] Yuhan Hao, Stephanie Hao, Erica Andersen-Nissen, William M Mauck, Shiwei Zheng, Andrew Butler, Maddie J Lee, Aaron J Wilk, Charlotte Darby, Michael Zager, et al.  Integrated analysis of multimodal single-cell data.  _Cell_ , 184(13):3573–3587, 2021.
  * He et al. [2016] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 770–778, 2016.
  * He et al. [2020] Kaiming He, Haoqi Fan, Yuxin Wu, Saining Xie, and Ross Girshick.  Momentum contrast for unsupervised visual representation learning.  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , pages 9729–9738, 2020.
  * He et al. [2022] Kaiming He, Xinlei Chen, Saining Xie, Yanghao Li, Piotr Dollár, and Ross Girshick.  Masked autoencoders are scalable vision learners.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 16000–16009, 2022.
  * Hou et al. [2022] Zhenyu Hou, Xiao Liu, Yukuo Cen, Yuxiao Dong, Hongxia Yang, Chunjie Wang, and Jie Tang.  Graphmae: Self-supervised masked graph autoencoders.  In _Proceedings of the 28th ACM SIGKDD Conference on Knowledge Discovery and Data Mining_ , pages 594–604, 2022.
  * Hu et al. [2020] Weihua Hu, Matthias Fey, Marinka Zitnik, Yuxiao Dong, Hongyu Ren, Bowen Liu, Michele Catasta, and Jure Leskovec.  Open graph benchmark: Datasets for machine learning on graphs.  _Advances in neural information processing systems_ , 33:22118–22133, 2020.
  * Huang et al. [2021] Yu Huang, Chenzhuang Du, Zihui Xue, Xuanyao Chen, Hang Zhao, and Longbo Huang.  What makes multi-modal learning better than single (provably).  _Advances in Neural Information Processing Systems_ , 34:10944–10956, 2021.
  * Hwang et al. [2018] Byungjin Hwang, Ji Hyun Lee, and Duhee Bang.  Single-cell rna sequencing technologies and bioinformatics pipelines.  _Experimental & molecular medicine_, 50(8):1–14, 2018.
  * Khromov and Singh [2023] Grigory Khromov and Sidak Pal Singh.  Some fundamental aspects about lipschitz continuity of neural network functions.  _arXiv preprint arXiv:2302.10886_ , 2023.
  * Kipf and Welling [2016] Thomas N Kipf and Max Welling.  Variational graph auto-encoders.  _NIPS Workshop on Bayesian Deep Learning_ , 2016.
  * Kipf and Welling [2017] Thomas N. Kipf and Max Welling.  Semi-supervised classification with graph convolutional networks.  In _International Conference on Learning Representations_ , 2017.  URL https://openreview.net/forum?id=SJU4ayYgl.
  * Krämer et al. [2014] Andreas Krämer, Jeff Green, Jack Pollard Jr, and Stuart Tugendreich.  Causal analysis approaches in ingenuity pathway analysis.  _Bioinformatics_ , 30(4):523–530, 2014.
  * Krassowski [2020] M Krassowski.  Complexupset: Create complex upset plots using ggplot2 components.  _R package version 0.5_ , 18, 2020.
  * Langfelder and Horvath [2008] Peter Langfelder and Steve Horvath.  Wgcna: an r package for weighted correlation network analysis.  _BMC bioinformatics_ , 9(1):1–13, 2008.
  * Lee et al. [2022] Patty J Lee, Philip Blood, Katy Börner, Judith Campisi, Feng Chen, Heike Daldrup-Link, Phil De Jager, Li Ding, Francesca E Duncan, Oliver Eickelberg, et al.  Nih sennet consortium: Mapping senescent cells in the human body to understand health and disease, 2022.
  * Leek et al. [2010] Jeffrey T Leek, Robert B Scharpf, Héctor Corrada Bravo, David Simcha, Benjamin Langmead, W Evan Johnson, Donald Geman, Keith Baggerly, and Rafael A Irizarry.  Tackling the widespread and critical impact of batch effects in high-throughput data.  _Nature Reviews Genetics_ , 11(10):733–739, 2010.
  * Liang et al. [2022] Paul Pu Liang, Amir Zadeh, and Louis-Philippe Morency.  Foundations and recent trends in multimodal machine learning: Principles, challenges, and open questions.  _arXiv preprint arXiv:2209.03430_ , 2022.
  * Lin et al. [2022] Xiang Lin, Tian Tian, Zhi Wei, and Hakon Hakonarson.  Clustering of single-cell multi-omics data with a multimodal deep learning method.  _Nature communications_ , 13(1):7705, 2022.
  * Liu et al. [2008] Yinyin Liu, Janusz A Starzyk, and Zhen Zhu.  Optimized approximation algorithm in neural networks without overfitting.  _IEEE transactions on neural networks_ , 19(6):983–995, 2008.
  * Lopez et al. [2018] Romain Lopez, Jeffrey Regier, Michael B Cole, Michael I Jordan, and Nir Yosef.  Deep generative modeling for single-cell transcriptomics.  _Nature methods_ , 15(12):1053–1058, 2018.
  * Luecken et al. [2022] Malte D Luecken, Maren Büttner, Kridsadakorn Chaichoompu, Anna Danese, Marta Interlandi, Michaela F Müller, Daniel C Strobl, Luke Zappia, Martin Dugas, Maria Colomé-Tatché, et al.  Benchmarking atlas-level data integration in single-cell genomics.  _Nature methods_ , 19(1):41–50, 2022.
  * Luo et al. [2017] Chongyuan Luo, Christopher L Keown, Laurie Kurihara, Jingtian Zhou, Yupeng He, Junhao Li, Rosa Castanon, Jacinta Lucero, Joseph R Nery, Justin P Sandoval, et al.  Single-cell methylomes identify neuronal subtypes and regulatory elements in mammalian cortex.  _Science_ , 357(6351):600–604, 2017.
  * Ma et al. [2021] Kaili Ma, Haochen Yang, Han Yang, Tatiana Jin, Pengfei Chen, Yongqiang Chen, Barakeel Fanseu Kamhoua, and James Cheng.  Improving graph representation learning by contrastive regularization.  _arXiv preprint arXiv:2101.11525_ , 2021.
  * Marguerat and Bähler [2010] Samuel Marguerat and Jürg Bähler.  Rna-seq: from technology to biology.  _Cellular and molecular life sciences_ , 67:569–579, 2010.
  * McInnes et al. [2018] Leland McInnes, John Healy, Nathaniel Saul, and Lukas Großberger.  Umap: Uniform manifold approximation and projection.  _Journal of Open Source Software_ , 3(29):861, 2018.  doi: 10.21105/joss.00861.  URL https://doi.org/10.21105/joss.00861.
  * Meral [2022] Bekir Fatih Meral.  Parental views of families of children with autism spectrum disorder and developmental disorders during the covid-19 pandemic.  _Journal of Autism and Developmental Disorders_ , 52(4):1712–1724, 2022.
  * Misra [2019] Diganta Misra.  Mish: A self regularized non-monotonic activation function.  _arXiv preprint arXiv:1908.08681_ , 2019.
  * Mo et al. [2022] Yujie Mo, Liang Peng, Jie Xu, Xiaoshuang Shi, and Xiaofeng Zhu.  Simple unsupervised graph representation learning.  _Proceedings of the AAAI Conference on Artificial Intelligence_ , 36(7):7797–7805, Jun. 2022.  doi: 10.1609/aaai.v36i7.20748.  URL https://ojs.aaai.org/index.php/AAAI/article/view/20748.
  * Mohri et al. [2018] Mehryar Mohri, Afshin Rostamizadeh, and Ameet Talwalkar.  _Foundations of machine learning_.  MIT press, 2018.
  * Mulqueen et al. [2018] Ryan M Mulqueen, Dmitry Pokholok, Steven J Norberg, Kristof A Torkenczy, Andrew J Fields, Duanchen Sun, John R Sinnamon, Jay Shendure, Cole Trapnell, Brian J O’Roak, et al.  Highly scalable generation of dna methylation profiles in single cells.  _Nature biotechnology_ , 36(5):428–431, 2018.
  * Needham et al. [2020] Edward J Needham, Sherry H-Y Chou, Alasdair J Coles, and David K Menon.  Neurological implications of covid-19 infections.  _Neurocritical care_ , 32:667–671, 2020.
  * Nikzad et al. [2019] Rana Nikzad, Laura S Angelo, Kevin Aviles-Padilla, Duy T Le, Vipul K Singh, Lynn Bimler, Milica Vukmanovic-Stejic, Elena Vendrame, Thanmayi Ranganath, Laura Simpson, et al.  Human natural killer cells mediate adaptive immunity to viral antigens.  _Science immunology_ , 4(35):eaat8116, 2019.
  * Oh et al. [2022] Sehyun Oh, Ludwig Geistlinger, Marcel Ramos, Daniel Blankenberg, Marius van den Beek, Jaclyn N Taroni, Vincent J Carey, Casey S Greene, Levi Waldron, and Sean Davis.  Genomicsupersignature facilitates interpretation of rna-seq experiments through robust, efficient comparison to public databases.  _Nature Communications_ , 13(1):3695, 2022.
  * Oord et al. [2018] Aaron van den Oord, Yazhe Li, and Oriol Vinyals.  Representation learning with contrastive predictive coding.  _arXiv preprint arXiv:1807.03748_ , 2018.
  * Palit et al. [2019] Subarna Palit, Christoph Heuser, Gustavo P De Almeida, Fabian J Theis, and Christina E Zielinski.  Meeting the challenges of high-dimensional single-cell data analysis in immunology.  _Frontiers in immunology_ , 10:1515, 2019.
  * Pedregosa et al. [2011] Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, et al.  Scikit-learn: Machine learning in python.  _the Journal of machine Learning research_ , 12:2825–2830, 2011.
  * Rampášek et al. [2022] Ladislav Rampášek, Michael Galkin, Vijay Prakash Dwivedi, Anh Tuan Luu, Guy Wolf, and Dominique Beaini.  Recipe for a general, powerful, scalable graph transformer.  _Advances in Neural Information Processing Systems_ , 35:14501–14515, 2022.
  * Rosen et al. [2023] Yanay Rosen, Maria Brbić, Yusuf Roohani, Kyle Swanson, Ziang Li, and Jure Leskovec.  Towards universal cell embeddings: Integrating single-cell rna-seq datasets across species with saturn.  _Biorxiv: the Preprint Server for Biology_ , 2023.
  * Rozenblatt-Rosen et al. [2017] Orit Rozenblatt-Rosen, Michael JT Stubbington, Aviv Regev, and Sarah A Teichmann.  The human cell atlas: from vision to reality.  _Nature_ , 550(7677):451–453, 2017.
  * Saliba et al. [2014] Antoine-Emmanuel Saliba, Alexander J Westermann, Stanislaw A Gorski, and Jörg Vogel.  Single-cell rna-seq: advances and future challenges.  _Nucleic acids research_ , 42(14):8845–8860, 2014.
  * Salih et al. [2023] Amanda Salih, Aaron Chin, Manisha Gandhi, Amir Shamshirsaz, Hennie Lombaard, and Joud Hajjar.  Hereditary angioedema and covid-19 during pregnancy: Two case reports.  _The Journal of Allergy and Clinical Immunology: In Practice_ , 11(3):961–962, 2023.
  * Sassoulas et al. [2018] Pierre Sassoulas, Anazalea, and Marcomanz.  pySankey, 2018.  URL https://github.com/anazalea/pySankey.
  * Shi et al. [2021] Yunsheng Shi, Zhengjie Huang, Shikun Feng, Hui Zhong, Wenjing Wang, and Yu Sun.  Masked label prediction: Unified message passing model for semi-supervised classification.  In Zhi-Hua Zhou, editor, _Proceedings of the Thirtieth International Joint Conference on Artificial Intelligence, IJCAI-21_ , pages 1548–1554. International Joint Conferences on Artificial Intelligence Organization, 8 2021.  doi: 10.24963/ijcai.2021/214.  URL https://doi.org/10.24963/ijcai.2021/214.  Main Track.
  * Skinnider et al. [2019] Michael A Skinnider, Jordan W Squair, and Leonard J Foster.  Evaluating measures of association for single-cell transcriptomics.  _Nature methods_ , 16(5):381–386, 2019.
  * Srivastava et al. [2020] Rajneesh Srivastava, Swapna Vidhur Daulatabad, Mansi Srivastava, and Sarath Chandra Janga.  Role of sars-cov-2 in altering the rna-binding protein and mirna-directed post-transcriptional regulatory networks in humans.  _International journal of molecular sciences_ , 21(19):7090, 2020.
  * Sriwijitalai and Wiwanitkit [2020] Won Sriwijitalai and Viroj Wiwanitkit.  Hearing loss and covid-19: a note.  _American journal of otolaryngology_ , 2020.
  * Stoeckius et al. [2017] Marlon Stoeckius, Christoph Hafemeister, William Stephenson, Brian Houck-Loomis, Pratip K Chattopadhyay, Harold Swerdlow, Rahul Satija, and Peter Smibert.  Simultaneous epitope and transcriptome measurement in single cells.  _Nature methods_ , 14(9):865–868, 2017.
  * Strober et al. [2019] BJ Strober, Reem Elorbany, K Rhodes, Nirmal Krishnan, Karl Tayeb, Alexis Battle, and Yoav Gilad.  Dynamic genetic regulation of gene expression during cellular differentiation.  _Science_ , 364(6447):1287–1290, 2019.
  * Stuart et al. [2019] Tim Stuart, Andrew Butler, Paul Hoffman, Christoph Hafemeister, Efthymia Papalexi, William M Mauck, Yuhan Hao, Marlon Stoeckius, Peter Smibert, and Rahul Satija.  Comprehensive integration of single-cell data.  _Cell_ , 177(7):1888–1902, 2019.
  * Su et al. [2023] Chang Su, Zichun Xu, Xinning Shan, Biao Cai, Hongyu Zhao, and Jingfei Zhang.  Cell-type-specific co-expression inference from single cell rna-sequencing data.  _Nature Communications_ , 14(1):4846, 2023.
  * Sugiyama [2007] Masashi Sugiyama.  Dimensionality reduction of multimodal labeled data by local fisher discriminant analysis.  _Journal of machine learning research_ , 8(5), 2007.
  * Theodoris et al. [2023] Christina V Theodoris, Ling Xiao, Anant Chopra, Mark D Chaffin, Zeina R Al Sayed, Matthew C Hill, Helene Mantineo, Elizabeth M Brydon, Zexian Zeng, X Shirley Liu, et al.  Transfer learning enables predictions in network biology.  _Nature_ , pages 1–9, 2023.
  * Traag et al. [2019] Vincent A Traag, Ludo Waltman, and Nees Jan Van Eck.  From louvain to leiden: guaranteeing well-connected communities.  _Scientific reports_ , 9(1):1–12, 2019.
  * Tripuraneni et al. [2020] Nilesh Tripuraneni, Michael Jordan, and Chi Jin.  On the theory of transfer learning: The importance of task diversity.  _Advances in neural information processing systems_ , 33:7852–7862, 2020.
  * Tripuraneni et al. [2021] Nilesh Tripuraneni, Chi Jin, and Michael Jordan.  Provable meta-learning of linear representations.  In _International Conference on Machine Learning_ , pages 10434–10443. PMLR, 2021.
  * Van Dijk et al. [2018] David Van Dijk, Roshan Sharma, Juozas Nainys, Kristina Yim, Pooja Kathail, Ambrose J Carr, Cassandra Burdziak, Kevin R Moon, Christine L Chaffer, Diwakar Pattabiraman, et al.  Recovering gene interactions from single-cell data using data diffusion.  _Cell_ , 174(3):716–729, 2018.
  * Vaswani et al. [2017] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  _Advances in neural information processing systems_ , 30, 2017.
  * Veličković et al. [2018] Petar Veličković, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Liò, and Yoshua Bengio.  Graph attention networks.  In _International Conference on Learning Representations_ , 2018.  URL https://openreview.net/forum?id=rJXMpikCZ.
  * Virtanen et al. [2020] Pauli Virtanen, Ralf Gommers, Travis E Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, et al.  Scipy 1.0: fundamental algorithms for scientific computing in python.  _Nature methods_ , 17(3):261–272, 2020.
  * Wang et al. [2021a] Guangtao Wang, Rex Ying, Jing Huang, and Jure Leskovec.  Multi-hop attention graph neural networks.  In Zhi-Hua Zhou, editor, _Proceedings of the Thirtieth International Joint Conference on Artificial Intelligence, IJCAI-21_ , pages 3089–3096. International Joint Conferences on Artificial Intelligence Organization, 8 2021a.  doi: 10.24963/ijcai.2021/425.  URL https://doi.org/10.24963/ijcai.2021/425.  Main Track.
  * Wang et al. [2021b] Juexin Wang, Anjun Ma, Yuzhou Chang, Jianting Gong, Yuexu Jiang, Ren Qi, Cankun Wang, Hongjun Fu, Qin Ma, and Dong Xu.  scgnn is a novel graph neural network framework for single-cell rna-seq analyses.  _Nature communications_ , 12(1):1–11, 2021b.
  * Wang et al. [2020] Yue Wang, Shafiq Joty, Michael Lyu, Irwin King, Caiming Xiong, and Steven C.H. Hoi.  VD-BERT: A Unified Vision and Dialog Transformer with BERT.  In _Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing (EMNLP)_ , pages 3325–3338, Online, November 2020. Association for Computational Linguistics.  doi: 10.18653/v1/2020.emnlp-main.269.  URL https://aclanthology.org/2020.emnlp-main.269.
  * Wang et al. [2022] Yuge Wang, Tianyu Liu, and Hongyu Zhao.  Respan: a powerful batch correction model for scrna-seq data through residual adversarial networks.  _Bioinformatics_ , 38(16):3942–3949, 2022.
  * Wickham [2011] Hadley Wickham.  ggplot2.  _Wiley interdisciplinary reviews: computational statistics_ , 3(2):180–185, 2011.
  * Woicik et al. [2023] Addie Woicik, Mingxin Zhang, Hanwen Xu, Sara Mostafavi, and Sheng Wang.  Gemini: memory-efficient integration of hundreds of gene networks with high-order pooling.  _Bioinformatics_ , 39, 06 2023.  doi: 10.1093/bioinformatics/btad247.
  * Wolf et al. [2018] F Alexander Wolf, Philipp Angerer, and Fabian J Theis.  Scanpy: large-scale single-cell gene expression data analysis.  _Genome biology_ , 19(1):1–5, 2018.
  * Xie et al. [2021] Renjie Xie, Wei Xu, Yanzhi Chen, Jiabao Yu, Aiqun Hu, Derrick Wing Kwan Ng, and A. Lee Swindlehurst.  A generalizable model-and-data driven approach for open-set rff authentication.  _IEEE Transactions on Information Forensics and Security_ , 16:4435–4450, 2021.  doi: 10.1109/TIFS.2021.3106166.
  * Xu et al. [2019] Keyulu Xu, Weihua Hu, Jure Leskovec, and Stefanie Jegelka.  How powerful are graph neural networks?  In _International Conference on Learning Representations_ , 2019.  URL https://openreview.net/forum?id=ryGs6iA5Km.
  * Xu et al. [2022] Mingcong Xu, Xuefeng Bai, Bo Ai, Guorui Zhang, Chao Song, Jun Zhao, Yuezhu Wang, Ling Wei, Fengcui Qian, Yanyu Li, et al.  Tf-marker: a comprehensive manually curated database for transcription factors and related markers in specific cell and tissue types in human.  _Nucleic acids research_ , 50(D1):D402–D412, 2022.
  * Yang et al. [2022] Fan Yang, Wenchuan Wang, Fang Wang, Yuan Fang, Duyu Tang, Junzhou Huang, Hui Lu, and Jianhua Yao.  scbert as a large-scale pretrained deep language model for cell type annotation of single-cell rna-seq data.  _Nature Machine Intelligence_ , 4(10):852–866, 2022.
  * Yang et al. [2021] Karren Dai Yang, Anastasiya Belyaeva, Saradha Venkatachalapathy, Karthik Damodaran, Abigail Katcoff, Adityanarayanan Radhakrishnan, GV Shivashankar, and Caroline Uhler.  Multi-domain translation between single-cell imaging and sequencing data using autoencoders.  _Nature communications_ , 12(1):1–10, 2021.
  * Ying et al. [2021] Chengxuan Ying, Tianle Cai, Shengjie Luo, Shuxin Zheng, Guolin Ke, Di He, Yanming Shen, and Tie-Yan Liu.  Do transformers really perform badly for graph representation?  _Advances in Neural Information Processing Systems_ , 34:28877–28888, 2021.
  * Zhang et al. [2023a] Bohang Zhang, Shengjie Luo, Liwei Wang, and Di He.  Rethinking the expressive power of GNNs via graph biconnectivity.  In _The Eleventh International Conference on Learning Representations_ , 2023a.  URL https://openreview.net/forum?id=r9hNv76KoT3.
  * Zhang et al. [2023b] Di Zhang, Yanxiang Deng, Petra Kukanja, Eneritz Agirre, Marek Bartosovic, Mingze Dong, Cong Ma, Sai Ma, Graham Su, Shuozhen Bao, et al.  Spatial epigenome–transcriptome co-profiling of mammalian tissues.  _Nature_ , pages 1–10, 2023b.
  * Zhang et al. [2023c] Shuang Zhang, Rui Fan, Yuti Liu, Shuang Chen, Qiao Liu, and Wanwen Zeng.  Applications of transformer-based language models in bioinformatics: A survey.  _Bioinformatics Advances_ , 2023c.
  * Zheng et al. [2017] Grace XY Zheng, Jessica M Terry, Phillip Belgrader, Paul Ryvkin, Zachary W Bent, Ryan Wilson, Solongo B Ziraldo, Tobias D Wheeler, Geoff P McDermott, Junjie Zhu, et al.  Massively parallel digital transcriptional profiling of single cells.  _Nature communications_ , 8(1):1–12, 2017.
  * Zhou et al. [2022] Kaixiong Zhou, Xiao Huang, Qingquan Song, Rui Chen, and Xia Hu.  Auto-gnn: Neural architecture search of graph neural networks.  _Frontiers in Big Data_ , 5, 2022.  ISSN 2624-909X.  doi: 10.3389/fdata.2022.1029307.  URL https://www.frontiersin.org/articles/10.3389/fdata.2022.1029307.
  * Zhu et al. [2021a] Jiaqiang Zhu, Shiquan Sun, and Xiang Zhou.  Spark-x: non-parametric modeling enables scalable and robust detection of spatial expression patterns for large spatial transcriptomic studies.  _Genome Biology_ , 22(1):1–25, 2021a.
  * Zhu et al. [2021b] Jiong Zhu, Ryan A. Rossi, Anup Rao, Tung Mai, Nedim Lipka, Nesreen K. Ahmed, and Danai Koutra.  Graph neural networks with heterophily.  _Proceedings of the AAAI Conference on Artificial Intelligence_ , 35(12):11168–11176, May 2021b.  doi: 10.1609/aaai.v35i12.17332.  URL https://ojs.aaai.org/index.php/AAAI/article/view/17332.
  * Zhu et al. [2020] Yanqiao Zhu, Yichen Xu, Feng Yu, Qiang Liu, Shu Wu, and Liang Wang.  Deep graph contrastive representation learning.  _arXiv preprint arXiv:2006.04131_ , 2020.
  * Zitnik and Leskovec [2017] Marinka Zitnik and Jure Leskovec.  Predicting multicellular function through multi-layer tissue networks.  _Bioinformatics_ , 33(14):i190–i198, 2017.

## Appendix

###  Appendix A GIANT’s selection of datasets

Figure 5: Problematic datasets selected by GIANT.

The quality of datasets, as depicted in Figure 5, impacts the conclusions
derived from gene embedding learning methods. GIANT employed low-quality
spatial transcriptomic and scATAC-seq data, which ideally should be
circumvented at the research inception. The spatial data, HBM477.BPMF.855 and
HBM983.ZLKS.328, exhibit considerable barcode expression measurement loss.
Under normal circumstances, gene reads distribution in sectioned samples
should form complete circles instead of truncated patterns. The scATAC-seq
datasets, HBM737.XZCH.767 and HBM793.TLNJ.249, present challenges in
discerning differences in gene expression assays across distinct cells, as
evidenced by the UMAP visualization results. In an ideal scenario, cells
marked with varying colors should occupy different positions in the low-
dimensional representation. HBM628.XTDP.423 and HBM292.XZZJ.243 serve as
examples of high-quality datasets.

###  Appendix B Constructing co-expression network and training datasets

In this section, we outline the algorithmic details of preprocessing steps and
network construction using CS-CORE, scTransform, and SPARK-X. Numerous
databases exist for atlas-scale transcriptomic and epigenomic data analysis,
including the Human BioMolecular Atlas Program (HuBMAP) [28], the Human Cell
Atlas (HCA) [77], the 20 Cellular Senescence Network (SenNet) [53], and more.
We gather Single-cell RNA sequencing (scRNA-seq) from HCA datasets, [14] and
[76], Single-cell sequencing Assay for Transposase-accessible Chromatin
(scATAC-seq) from HuBMAP datasets and [14], and high-quality spatial data from
[3, 112]. For distinct omics data, we implement one common and one specific
process.

In the common process, we filter barcodes with gene expression counts below
200 and genes with expression counts below three in each barcode. We also
filter Mitochondrial (MT) genes. For scRNA-seq datasets, we employ scTransform
[35] to select HVGs and generate Pearson residuals, replacing raw expression
with residuals. scTransform is the first model to incorporate sequencing depth
as a covariate rather than directly applying the size factor to normalization.
It eliminates confounding caused by sequencing depth in raw single-cell or
spatial expression data, generating corrected gene expression profiles. These
advantages make it a widely used normalization method [11].

To construct scRNA-seq dataset co-expression networks based on Unique
Molecular Identifier (UMI), we use CS-CORE [88], a state-of-the-art tool for
co-expression inference based on UMI count data. CS-CORE demonstrates
increased robustness and a lower false positive rate compared to other tools.
For scATAC-seq datasets, we use Seurat [39] to convert the original cells-
peaks matrix into the cell-gene activity matrix, incorporating prior
information. The cell-gene activity matrix can be processed similarly to the
scRNA-seq data matrix, so subsequent preprocessing steps remain the same.

To construct spatial data co-expression networks, we consider spatial
expression patterns (SE genes or spatially HVGs) and treat each barcode as a
sample. We identify SE genes using SPARK-X [116], then generate corrected gene
expression profiles based on scTransform and co-expression networks based on
CS-CORE.

####  B.1 CS-CORE

Since our data are in UMI type, CS-CORE can be used to perform the co-
expression network construction. Now considering we have $n$ cells and for one
cell $i$, its expression profile can be denoted as a vector
$(x_{i1},...,x_{ip})$, where $p$ is the number of genes. We can also use
$s_{i}$ to represent the sequencing depth of cell $i$. Given the underlying
expression levels from cell $i$ with $p$ genes as
$\left(z_{i1},\ldots,z_{ip}\right)$, the assumption of CS-CORE for the
expression profile follows:

$\left(z_{i1},\ldots,z_{ip}\right)\sim
F_{p}(\boldsymbol{\mu},\boldsymbol{\Sigma}),\quad x_{ij}\mid
z_{ij}\sim\operatorname{Poisson}\left(s_{i}z_{ij}\right)$ (8)

where $F_{p}$ is an unknown non-negative p-variate distribution with
$\boldsymbol{\mu}$ as mean and $\boldsymbol{\Sigma}$ as covariance matrix.
Here the observed expression $x_{ij}$ follows the Poisson measurement model
depending on the underlying expression $z_{ij}$ and the sequencing depth
$s_{i}=\sum_{j=1}^{n}x_{ij}$. CS-CORE applies a moment-based iteratively
reweighted least squares (IRLS) estimation procedure to estimate the
covariance matrix. Once we have the covariance matrix
$\boldsymbol{\Sigma}_{p\times p}=[\sigma_{ij}]_{i=1...p}^{j=1...p}$, we can
use the correlation
$\rho_{jj^{\prime}}=\frac{\sigma_{jj^{\prime}}}{\sqrt{\sigma_{jj}\sigma_{j^{\prime}j^{\prime}}}}$
to estimate the co-expression relation between gene $j$ and gene $j^{\prime}$.

With such an assumption, CS-CORE can model the relationship between gene $j$
and gene $j^{\prime}$ based on a statistical test:

$T_{jj^{\prime}}=\frac{\sum_{i}s_{i}^{2}\left(x_{ij}-s_{i}\mu_{j}\right)\left(x_{ij^{\prime}}-s_{i}\mu_{j^{\prime}}\right)g_{ijj^{\prime}}}{\sqrt{\sum_{i}s_{i}^{4}\left(s_{i}\mu_{j}+s_{i}^{2}\sigma_{jj}\right)\left(s_{i}\mu_{j^{\prime}}+s_{i}^{2}\sigma_{j^{\prime}j^{\prime}}\right)g_{ijj^{\prime}}^{2}}}$
(9)

The statistic $T_{jj^{\prime}}$ under the null hypothesis assume gene $j$ and
gene $j^{\prime}$ are independent, which means there is no edge between these
two genes. We chose to use p-value ($p<0.005$) to construct the edges in the
co-expression network.

####  B.2 scTransform

To remove the confounding effect of sequencing depth from the expression
level, we first process the single-cell data using scTransform to get the
Pearson residuals and then use the Pearson residuals as the initial embeddings
of different genes. By assuming that the UMI count data follow the negative
binomial distribution, for a given gene $g$ in cell $c$, we have:

$\displaystyle x_{gc}$
$\displaystyle\sim\mathrm{NB}\left(\mu_{gc},\theta_{g}\right)$ (10)
$\displaystyle\ln\mu_{gc}$ $\displaystyle=\beta_{g0}+\ln s_{c},$

scTransform regularizes $\theta$ as a function of gene means $\mu$, by
utilizing the Generalized Linear Model (GLM) with a log link function provided
by the above equation. Furthermore, we can estimate the unknown parameters and
calculate the Pearson residuals $Z_{gc}$ based on:

$\displaystyle Z_{gc}=\frac{x_{gc}-\mu_{gc}}{\sigma_{gc}}$ (11)
$\displaystyle\mu_{gc}=\exp(\beta_{g0}+\ln s_{c})$
$\displaystyle\sigma_{gc}=\sqrt{\mu_{gc}+\frac{\mu_{gc}^{2}}{\theta_{gc}}}$

We can finally replace the original expression matrix with the residual matrix
$Z$ generated by GLM, and store the expression matrix and the corresponding
graph in Scanpy files.

####  B.3 SPARK-X

The primary distinctions between single-cell transcriptomic data and spatial
transcriptomic data involve two aspects: 1. In most spatially resolved data,
each barcode represents a mixture of different cells. 2. The additional
spatial information introduces spatial gene expression patterns (SE genes) for
spatially resolved data. To identify SE genes in spatial transcriptomic data,
SPARK-X employs a statistical test comparing the distance covariance matrix,
which is constructed using barcode positions, and the expression covariance
matrix, which is built from the gene expression profiles.

More specifically, for a spatial transcriptomic gene expression matrix with
size $n\times d$, we can donate the matrix for coordinates of samples as
$S=(s_{1}^{T},...,s_{n}^{T}),s_{i}=(s_{i1},s_{i2})$. Therefore, the whole
expression matrix can be represented as :
$y=(y_{1}(s_{1}),...,y_{n}(s_{n}))^{T}$. Our target is to test whether $y$ is
independent from $S$, so we construct the expression covariance matrix based
on $E=y(y^{T}y)^{-1}y^{T}$, and the distance covariance matrix based on
$\Sigma=S(S^{T}S)^{-1}S^{T}$. For the distance covariance matrix $S$, SPARK-X
considers different kernels to describe different spatial expression patterns,
including 1. Gaussian kernel
$(s^{\prime}_{i1},s^{\prime}_{i2})=(exp(\frac{-s^{2}_{i1}}{2\sigma_{1}^{2}}),exp(\frac{-s^{2}_{i2}}{2\sigma_{2}^{2}}))$
and 2. Cosine kernel $(s^{\prime}_{i1},s^{\prime}_{i2})=(cos(\frac{2\pi
s_{i1}}{\Phi_{1}}),cos(\frac{2\pi s_{i2}}{\Phi_{2}}))$. After centering these
two matrices, we have $E_{C}$ and $\Sigma_{C}$. Therefore, the test statistic
is:

$T=\frac{\operatorname{trace}(E_{C}\Sigma_{C})}{n}$ (12)

This statistic follows a $\chi^{2}_{1}$ distribution and we record the p-value
for each gene. The null hypothesis is that the gene expression is irrelevant
to the position of barcodes. After finding the p-value list, we rank the
p-value in ascending and select the top 1000 genes to perform normalization
and co-expression graph construction.

The whole process of the graph construction is shown in Figure 6.

Figure 6: Schematic diagram of the graph data construction methods.

####  B.4 Transformconv

If we consider $c$ as the index of attention heads, the multi-head attention
for GNN from node $j$ to node $i$ can be expressed as:

$\displaystyle q_{c,i}^{(l)}$
$\displaystyle=W_{c,q}^{(l)}h_{i}^{(l)}+b_{c,q}^{(l)};k_{c,j}^{(l)}=W_{c,k}^{(l)}h_{j}^{(l)}+b_{c,k}^{(l)},$
(13) $\displaystyle e_{c,ij}$
$\displaystyle=W_{c,e}e_{ij}+b_{c,e};\alpha_{c,ij}^{(l)}=\frac{\left\langle
q_{c,i}^{(l)},k_{c,j}^{(l)}+e_{c,ij}\right\rangle}{\sum_{u\in\mathcal{N}(i)}\left\langle
q_{c,i}^{(l)},k_{c,u}^{(l)}+e_{c,iu}\right\rangle},$

where $<q,k>=e^{\frac{q^{T}k}{\sqrt{d}}}$. $\sqrt{d}$ is a scalar used to
reduce gradient vanishing [95] as introduced in the original Transformer
paper. The different vectors $q,k$ correspond to the query vector and the key
vector, while $e$ represents the edge features. The attention
$\alpha^{(l)}_{c,ij}$ denotes the $c_{th}$ attention value from node $j$ to
$i$ for layer $l$. $h$ represents the node embedding,
$W_{c,q}^{(l)},W_{c,k}^{(l)},W_{c,e}$ are weight matrices, and
$b_{c,q}^{(l)},b_{c,k}^{(l)},b_{c,e}$ are bias terms.

We define the embedding of node $i$ in layer $l+1$ as $h_{i}^{l+1}$ and update
the node embedding by:

$\displaystyle v_{c,j}^{(l)}$
$\displaystyle=W_{c,v}^{(l)}h_{j}^{(l)}+b_{c,v}^{(l)},$ (14) $\displaystyle
h_{i}^{(l+1)}$
$\displaystyle=\|_{c=1}^{C}\left[\sum_{j\in\mathcal{N}(i)}\alpha_{c,ij}^{(l)}\left(v_{c,j}^{(l)}+e_{c,ij}\right)\right],$

where the vector $v$ represents the value vector. The operation $\|_{c=1}^{C}$
denotes the concatenation of $C$ total attention heads. $\mathcal{N}(i)$
represents the neighbors of node $i$. Additionally, we construct the residual
link [40] based on $h_{i}^{l-1}$ for $h_{i}^{l+1}$. The activation function
employed to connect different layers is Mish
($\text{Mish}(x)=\tanh(\ln(1+\exp(x)))$) [65]. We also incorporated a
GraphNorm layer [13] to the connection path between hidden layers to enhance
convergence.

###  Appendix C Differential co-expression analysis

In this section, we prove that using CS-CORE, we could obtain the co-
expression relationships from multimodal biological data without the negative
impact of confounding factors brought by the batch effect.

Theorem C.1 Given the true observed gene expression level $x$ and gene
expression affected by batch effect $x^{\prime}$, the construction of the co-
expression network is the same based on the statistical test.

Proof. We consider the point $x^{\prime}_{ij}=ax_{ij}+b$, where $a,b$
represent the effect caused by batch effect, and $x_{ij}$ represents the true
biological data. Moreover, based on the probabilistic model provided by CS-
CORE, we have
$\mathbb{E}\left(x_{ij}\right)=s_{i}\mu_{j},Var(x_{ij})=s_{i}\mu_{j}+s_{i}^{2}\sigma_{jj}$.
Therefore, for the two points $x_{ij}$ and $x^{\prime}_{ij}$ we have:

$\displaystyle x_{ij}=s_{i}\mu_{j}+\epsilon_{ij}$ (15)
$\displaystyle\left(x_{ij}-s_{i}\mu_{j}\right)^{2}=s_{i}\mu_{j}+s_{i}^{2}\sigma_{jj}+\eta_{ij}$
$\displaystyle\left(x_{ij}-s_{i}\mu_{j}\right)\left(x_{ij^{\prime}}-s_{i}\mu_{j^{\prime}}\right)=s_{i}^{2}\sigma_{jj^{\prime}}+\xi_{ijj^{\prime}}$

and

$\displaystyle x^{\prime}_{ij}=as_{i}\mu_{j}+b+a\epsilon_{ij}$ (16)
$\displaystyle\left(x^{\prime}_{ij}-(as_{i}\mu_{j}+b)\right)^{2}=a^{2}s_{i}\mu_{i}+s_{i}^{2}a^{2}\sigma_{jj}+\eta_{ij}$
$\displaystyle\left(x^{\prime}_{ij}-(as_{i}\mu_{j}+b)\right)\left(x^{\prime}_{ij^{\prime}}-(as_{i}\mu_{j^{\prime}}+b)\right)=s_{i}^{2}a^{2}\sigma_{jj^{\prime}}+\xi_{ijj^{\prime}}$

Based on this assumption, we can calculate the test statistic
$T^{\prime}_{jj^{\prime}}$ with batch effect, which is:

$\displaystyle T^{\prime}_{jj^{\prime}}$
$\displaystyle=\frac{\sum_{i}s_{i}^{2}\left(x^{\prime}_{ij}-as_{i}\mu_{j}-b\right)\left(x^{\prime}_{ij^{\prime}}-as_{i}\mu_{j^{\prime}}-b\right)g_{ijj^{\prime}}}{\sqrt{\sum_{i}s_{i}^{4}\left(a^{2}s_{i}\mu_{j}+a^{2}s_{i}^{2}\sigma_{jj}\right)\left(a^{2}s_{i}\mu_{j^{\prime}}+a^{2}s_{i}^{2}\sigma_{j^{\prime}j^{\prime}}\right)g_{ijj^{\prime}}^{2}}}$
(17)
$\displaystyle=\frac{a^{2}\sum_{i}s_{i}^{2}\left(x_{ij}-s_{i}\mu_{j}\right)\left(x_{ij^{\prime}}-s_{i}\mu_{j^{\prime}}\right)g_{ijj^{\prime}}}{a^{2}\sqrt{\sum_{i}s_{i}^{4}\left(s_{i}\mu_{j}+s_{i}^{2}\sigma_{jj}\right)\left(s_{i}\mu_{j^{\prime}}+s_{i}^{2}\sigma_{j^{\prime}j^{\prime}}\right)g_{ijj^{\prime}}^{2}}}$
$\displaystyle=\frac{\sum_{i}s_{i}^{2}\left(x_{ij}-s_{i}\mu_{j}\right)\left(x_{ij^{\prime}}-s_{i}\mu_{j^{\prime}}\right)g_{ijj^{\prime}}}{\sqrt{\sum_{i}s_{i}^{4}\left(s_{i}\mu_{j}+s_{i}^{2}\sigma_{jj}\right)\left(s_{i}\mu_{j^{\prime}}+s_{i}^{2}\sigma_{j^{\prime}j^{\prime}}\right)g_{ijj^{\prime}}^{2}}}$
$\displaystyle=T_{jj^{\prime}}$

This statistic is irrelevant to $a,b$. Therefore, the batch effect cannot
affect our co-expression calculation, and our correlation defined by CS-CORE
can reflect the correlation for true biological content.

###  Appendix D Theory analysis for multimodal machine learning

Here we present a series of explanation for why using more modalities is more
suitable in our task based on arguing that MMML method can generate more
accurate estimate of the latent space representation for genes as the number
of modalities increases.

Here we consider our datasets set $\mathcal{D}$ and different types of
modalities $\mathcal{M}$ and $\mathcal{N}$, and their corresponding gene
embeddings as $e_{\mathcal{M}}$ and $e_{\mathcal{N}}$, where
$\textsc{Card}(\mathcal{N})<\textsc{Card}(\mathcal{M})$. We also have the true
representation $e^{*}$. Based on our model, we can estimate
$\hat{e}_{\mathcal{M}},\hat{e}_{\mathcal{N}}$. Here we intend to prove the
theorem:

Theorem D.1 $\hat{e}_{\mathcal{M}}$ is a better estimation for $e^{*}$
comparing to $\hat{e}_{\mathcal{N}}$ with
$\textsc{Card}(\mathcal{N})<\textsc{Card}(\mathcal{M})$.

Here we consider dataset $\mathcal{S}=\\{(X_{i},Y_{i})\\}_{i=1}^{m}$ as a
general training dataset. In our specific cases, we can treat $X_{i}$ carries
information with graph structure and $Y_{i}$ contains the edges information of
$X_{i}$. Therefore, we have formalized the dataset we used in our unsupervised
learning task. Our target is to minimize the empirical risk:

$\begin{aligned} \min&\quad\hat{r}\left(h\circ
g_{\mathcal{M}}\right)\triangleq\frac{1}{m}\sum_{i=1}^{m}\ell\left(h\circ
g_{\mathcal{M}}\left(\mathbf{X}_{i}\right),Y_{i}\right)\\\ \text{ s.t. }&\quad
h\in\mathcal{H},g_{\mathcal{M}}\in\mathcal{G}_{\mathcal{M}}\end{aligned},$
(18)

where $l(\cdot,\cdot)$ represents the loss function, $h\circ
g_{\mathcal{M}}\left(\mathbf{X}_{i}\right)=h(g_{\mathcal{M}}(\mathbf{X}_{i}))$
represents the composite function of $h$ and $g_{\mathcal{M}}$. Here function
$g$ maps data from $X_{i}$ to the latent space, while function $h$ maps data
from latent space to the space of $Y_{i}$. We have $\mathcal{H}$ as the
function class of $h$ and $G_{\mathcal{M}}$ is the function class of
$g_{\mathcal{M}}$. We define the true map functions as $h^{*}$ and $g^{*}$. We
donate $r(\cdot,\cdot)$ as risk.

To finish our proof, we need three assumptions.

Assumption D.1. The loss function is L-smooth with respect to the first
ordinate, and is bounded by a constant C [67, 92, 93].

In our case, we have three components. BCELoss and Cosine Similarity follow
this property [105, 47]. InfoNCE loss also follows this property proved by
[61].

Assumption D.2. The true latent representation $g^{\star}$ is contained in
$\mathcal{G}$, and the task mapping $h^{\star}$ is contained in $\mathcal{H}$
[25, 92, 93].

Assumption D.3. For any $g^{\prime}\in\mathcal{G}^{\prime}$ and
$\mathcal{M}\subset[K],g^{\prime}\circ
p_{\mathcal{M}}^{\prime}\in\mathcal{G}^{\prime}$ [45].

Here the set $[K]$ represents modality set with $K$ modalities.
$p_{\mathcal{M}}^{\prime}$ is a diagonal matrix using 1 for $ii$-th entry
$\in\mathcal{M}$ and 0 otherwise.

We also need one metric to evaluate the quality of our generated embeddings,
so here we introduce our definition of latent representation quality.

Definition D.1. For any latent space generated by $g$, the latent
representation quality is defined as [45]

$\eta(g)=\inf_{h\in\mathcal{H}}\left[r(h\circ g)-r\left(h^{*}\circ
g^{*}\right)\right]$ (19)

One approach we can use to prove Theorem D.1 is to explore the order
relationship between the upper bound of $\eta(g_{\mathcal{M}})$ and the upper
bound of $\eta(g_{\mathcal{N}})$. That is, we intend to prove:

$\sup_{g_{\mathcal{M}}\in
G_{\mathcal{M}}}\eta(g_{\mathcal{M}})\leq\sup_{g_{\mathcal{N}}\in
G_{\mathcal{N}}}\eta(g_{\mathcal{N}})$

To prove the above relation, we rely on a new theorem to compute the upper
bound of $\eta(\cdot)$.

Theorem D.2. [45] Let
$\mathcal{S}=\left\\{\left(\mathbf{X}_{i},Y_{i}\right)\right\\}_{i=1}^{m}$ be
a dataset of m examples drawn i.i.d. according to $\mathcal{D}$. Let
$\mathcal{M}$ be a subset of $[K]$. Assuming we have produced the empirical
risk minimizers $\left(\hat{h}_{\mathcal{M}},\hat{g}_{\mathcal{M}}\right)$
training with the $\mathcal{M}$ modalities. Then, for all $1>\delta>0$ with
probability at least $1-\delta:$

$\eta\left(\hat{g}_{\mathcal{M}}\right)\leq
4L\Re_{m}\left(\mathcal{H}\circ\mathcal{G}_{\mathcal{M}}\right)+4L\Re_{m}(\mathcal{H}\circ\mathcal{G})+6C\sqrt{\frac{2\ln(2/\delta)}{m}}+\hat{L}\left(\hat{h}_{\mathcal{M}}\circ\hat{g}_{\mathcal{M}},\mathcal{S}\right),$

where
$\hat{L}\left(\hat{h}_{\mathcal{M}}\circ\hat{g}_{\mathcal{M}},\mathcal{S}\right)\triangleq\hat{r}\left(\hat{h}_{\mathcal{M}}\circ\hat{g}_{\mathcal{M}}\right)-\hat{r}\left(h^{\star}\circ
g^{\star}\right)$ is the centered empirical loss and $L$ represents the
smoothness coefficient.

Remark. Now we consider sets $\mathcal{N}\subset\mathcal{M}\subset[K]$, and
based on Assumption D.3, we have
$\mathcal{G}_{\mathcal{N}}\subset\mathcal{G}_{\mathcal{M}}\subset\mathcal{G}$.
We can interpret this relation as the size of parameter space Param follows
the relation:
$\textsc{Param}_{\mathcal{N}}\subset\textsc{Param}_{\mathcal{M}}\subset\textsc{Param}_{\mathcal{G}}$.
Therefore, larger function class has a smaller empirical risk, so we have

$\hat{L}\left(\hat{h}_{\mathcal{M}}\circ\hat{g}_{\mathcal{M}},\mathcal{S}\right)\leq\hat{L}\left(\hat{h}_{\mathcal{N}}\circ\hat{g}_{\mathcal{N}},\mathcal{S}\right)$
(20)

Based on Theorem D.2, we can derive the upper bound of
$\eta(\hat{g}_{\mathcal{N}})$ as:

$\eta\left(\hat{g}_{\mathcal{N}}\right)\leq
4L\Re_{m}\left(\mathcal{H}\circ\mathcal{G}_{\mathcal{N}}\right)+4L\Re_{m}(\mathcal{H}\circ\mathcal{G})+6C\sqrt{\frac{2\ln(2/\delta)}{m}}+\hat{L}\left(\hat{h}_{\mathcal{N}}\circ\hat{g}_{\mathcal{N}},\mathcal{S}\right),$

where $\Re_{m}(\cdot)$ represents the Rademacher complexity [10]. Here we can
use the Right-hand Side (RHS) to approximate term
$4L\Re_{m}\left(\mathcal{H}\circ\mathcal{G}_{\mathcal{M}}\right)\sim\sqrt{C\left(\mathcal{H}\circ\mathcal{G}_{\mathcal{M}}\right)/m}$
and term
$4L\Re_{m}\left(\mathcal{H}\circ\mathcal{G}_{\mathcal{N}}\right)\sim\sqrt{C\left(\mathcal{H}\circ\mathcal{G}_{\mathcal{N}}\right)/m}$.
Based on the basic structural property of Rademacher complexity [10], we have

$C\left(\mathcal{H}\circ\mathcal{G}_{\mathcal{N}}\right)\leq
C\left(\mathcal{H}\circ\mathcal{G}_{\mathcal{M}}\right)$

Therefore, we can compute:

$\displaystyle\eta\left(\hat{g}_{\mathcal{M}}\right)-\eta\left(\hat{g}_{\mathcal{N}}\right)$
$\displaystyle=4L\Re_{m}\left(\mathcal{H}\circ\mathcal{G}_{\mathcal{M}}\right)+\hat{L}\left(\hat{h}_{\mathcal{M}}\circ\hat{g}_{\mathcal{M}},\mathcal{S}\right)-4L\Re_{m}\left(\mathcal{H}\circ\mathcal{G}_{\mathcal{N}}\right)-\hat{L}\left(\hat{h}_{\mathcal{N}}\circ\hat{g}_{\mathcal{N}},\mathcal{S}\right)$
$\displaystyle=\sqrt{\frac{C\left(\mathcal{H}\circ\mathcal{G}_{\mathcal{M}}\right)}{m}}-\sqrt{\frac{C\left(\mathcal{H}\circ\mathcal{G}_{\mathcal{N}}\right)}{m}}+\hat{L}\left(\hat{h}_{\mathcal{M}}\circ\hat{g}_{\mathcal{M}},\mathcal{S}\right)-\hat{L}\left(\hat{h}_{\mathcal{N}}\circ\hat{g}_{\mathcal{N}},\mathcal{S}\right)$
$\displaystyle=\frac{\sqrt{C\left(\mathcal{H}\circ\mathcal{G}_{\mathcal{M}}\right)}-\sqrt{C\left(\mathcal{H}\circ\mathcal{G}_{\mathcal{N}}\right)}}{\sqrt{m}}+\hat{L}\left(\hat{h}_{\mathcal{M}}\circ\hat{g}_{\mathcal{M}},\mathcal{S}\right)-\hat{L}\left(\hat{h}_{\mathcal{N}}\circ\hat{g}_{\mathcal{N}},\mathcal{S}\right)$
$\displaystyle\rightarrow\hat{L}\left(\hat{h}_{\mathcal{M}}\circ\hat{g}_{\mathcal{M}},\mathcal{S}\right)-\hat{L}\left(\hat{h}_{\mathcal{N}}\circ\hat{g}_{\mathcal{N}},\mathcal{S}\right)|_{m\rightarrow\infty}$
$\displaystyle\leq 0$

Remark. Therefore, as the number of modalities increases, the representation
quality of more modalities will be better than the representation quality of
fewer modalities. Hence we proved Theorem D.1 under large number of modalities
cases. In the real biological application, genes can be active in many
different scenarios, including different cells, different tissues and even
different species. Therefore, our proof fits the context of MuSe-GNN’s real
application.

###  Appendix E Extra experiments, ablation test and model selection

####  E.1 Ablation tests

We further investigate the contributions of various loss functions and GNN
models to MuSe-GNN’s performance. Our final loss function comprises three
components: the graph reconstruction loss $\mathcal{L}_{\text{reconst}}$, the
similarity maximization loss $\mathcal{L}_{\text{CosSim}}$, and the
contrastive learning loss $\mathcal{L}_{\text{InfoNCE}}$. The reconstruction
loss serves as the basic loss, and we assess our model’s performance with
different loss function components across all scRNA-seq datasets. The results
are presented in Table 3 and 4.

Table 3: The overall benchmarking average rank for ablation tests based on different tissues. $\mathcal{L}_{\text{CosSim}}$ | $\mathcal{L}_{\text{InfoNCE}}$ | Heart | Lung | Liver | Kidney | Thymus | Spleen | Pancreas | Cerebrum | Cerebellum | PBMC | Avg Rank
---|---|---|---|---|---|---|---|---|---|---|---|---
|  | 2.83 | 3.50 | 3.17 | 3.33 | 2.67 | 3.17 | 3.00 | 3.67 | 3.00 | 3.33 | 3.17
✓ |  | 2.17 | 1.67 | 2.00 | 2.17 | 2.50 | 2.00 | 1.50 | 2.00 | 2.17 | 5.67 | 2.38
| ✓ | 3.00 | 3.00 | 3.17 | 2.83 | 3.17 | 3.17 | 3.33 | 3.00 | 2.83 | 5.00 | 3.25
✓ | ✓ | 2.00 | 1.83 | 1.67 | 1.67 | 1.67 | 1.67 | 2.17 | 1.33 | 2.00 | 3.83 | 1.98
Table 4: The overall benchmarking average score for ablation tests based on different tissues. $\mathcal{L}_{\text{CosSim}}$ | $\mathcal{L}_{\text{InfoNCE}}$ | Heart | Lung | Liver | Kidney | Thymus | Spleen | Pancreas | Cerebrum | Cerebellum | PBMC | Avg Score
---|---|---|---|---|---|---|---|---|---|---|---|---
|  | 0.42 | 0.13 | 0.21 | 0.24 | 0.36 | 0.19 | 0.28 | 0.09 | 0.37 | 0.23 | 0.25
✓ |  | 0.65 | 0.81 | 0.77 | 0.75 | 0.57 | 0.70 | 0.83 | 0.71 | 0.68 | 0.66 | 0.71
| ✓ | 0.34 | 0.27 | 0.25 | 0.23 | 0.23 | 0.20 | 0.18 | 0.26 | 0.35 | 0.23 | 0.25
✓ | ✓ | 0.78 | 0.81 | 0.94 | 0.80 | 0.70 | 0.96 | 0.53 | 0.91 | 0.74 | 0.90 | 0.81

Based on the results of this ablation test, we conclude that combining all
three components to construct our loss function is the optimal choice,
achieving a performance that is 224.0% higher than the version without any
additional regularization. Furthermore, using only contrastive learning may
reduce the performance of our model, implying that learning similarity is more
crucial than learning differences for the gene embedding generation task. When
both components are included, the final average rank is 14.1% higher than the
version with only the similarity learning part. Thus, we need to incorporate
both weighted similarity learning and contrastive learning in our final design
to learn unified gene embeddings with improved performance.

Additionally, we compare the TransformConv framework with other GNN models,
including GCN, Graph Attention Network (GAT) [96], SUGRL (GCN+MLP+Contrastive
Learning) [66], GPS [75], GRACE (GCN+Contrastive Learning) [118], GraphSAGE
[36], Graph Isomorphism Network (GIN) [106], GraphMAE [43] and Graphormer
[110]. We set the number of epochs to 1000 for this comparison. The results
are displayed in Table 5.

Table 5: The overall benchmarking score and rank for model selection. Methods | ASW | AUC | iLISI | GC | CGR | NO | Avg Rank | Avg Score
---|---|---|---|---|---|---|---|---
GCN |  0.74$\pm$0.02 |  0.82$\pm$0.01 |  0.32$\pm$0.02 |  0.44$\pm$0.02 |  0.32$\pm$0.04 |  0.31$\pm$0.01 | 3.5 | 0.51
GAT |  0.74$\pm$0.09 | 0.83$\pm$0.01 |  0.14$\pm$0.03 |  0.42$\pm$0.03 |  0.00$\pm$0.00 |  0.15$\pm$0.03 | 5 | 0.17
SUGRL | 0.88$\pm$0.01 |  0.62$\pm$0.00 |  0.45$\pm$0.01 |  0.43$\pm$0.03 |  0.38$\pm$0.02 |  0.29$\pm$0.00 | 3.5 | 0.56
GPS |  0.77$\pm$0.02 |  0.63$\pm$0.00 |  0.50$\pm$0.01 | 0.72$\pm$0.03 |  0.57$\pm$0.02 |  0.30$\pm$0.01 | 2.67 | 0.68
GRACE |  0.82$\pm$0.02 |  0.81$\pm$0.00 |  0.35$\pm$0.01 |  0.42$\pm$0.04 |  0.11$\pm$0.01 |  0.25$\pm$0.01 | 4.17 | 0.48
GraphSAGE | OOM | OOM | OOM | OOM | OOM | OOM | 7.00 | 0
GIN | OOM | OOM | OOM | OOM | OOM | OOM | 7.00 | 0
GraphMAE | OOM | OOM | OOM | OOM | OOM | OOM | 7.00 | 0
Graphormer | OOM | OOM | OOM | OOM | OOM | OOM | 7.00 | 0
Transformconv |  0.75$\pm$0.01 |  0.78$\pm$0.02 | 0.51$\pm$0.02 |  0.68$\pm$0.04 | 0.61$\pm$0.02 | 0.31$\pm$0.00 | 2.17 | 0.80

In this table, OOM means out of memory. We can conclude that using
TransformConv as the basic graph neural network framework is the best choice,
which is 17.6% higher than the version based on GPS (the top2 model).
Moreover, based on Figure 7, we can also discover that methods based on GCN or
GAT failed to learn the gene function similarity across different datasets.
Therefore, choosing TransformConv is reasonable, and the major contribution of
Transformer model towards gene embedding learning task is the multi-head
attention design.

Figure 7: Gene embeddings based on different GNNs.

We also conduct experiments to justify the need for weight-sharing design. In
this comparison, we evaluate the performance and model statistics of two
models (with weight-sharing design (ws) and without weight-sharing design (w/o
ws)). As shown in Table 6, the combined high or equivalent rank of weight-
sharing designs in comparison to w/o weight-sharing designs is nearly
consistent across all tissue types. Furthermore, the weight-sharing design can
also reduce the model’s training cost by decreasing the model size, the number
of trainable parameters, and the total training time, as confirmed in Table 7.

Table 6: Comparison between with weight-sharing design and without weight-sharing design. Methods | Heart | Lung | Liver | Kidney | Thymus | Spleen | Pancreas | Cerebrum | Cerebellum | PBMC | Avg Rank
---|---|---|---|---|---|---|---|---|---|---|---
ws | 1.50 | 1.33 | 1.33 | 1.33 | 1.33 | 1.33 | 1.67 | 1.50 | 1.17 | 1.17 | 1.37
w/o ws | 1.50 | 1.67 | 1.67 | 1.67 | 1.67 | 1.67 | 1.33 | 1.50 | 1.83 | 1.83 | 1.63
Table 7: Model statistics. Methods | model size (MB) | # of parameters (M) | training time per epoch (s)
---|---|---|---
ws | 1333 | 349 | 3.64
w/o ws | 1346.8 | 353 | 4.85

To demonstrate the necessity of both expression profiles and graph structure
for learning gene embeddings with biological information, we design an
experiment in which all nodes in different graphs are replaced with the same
features while maintaining the original graph structure. Based on the results
shown in Figure 8 and Table 8, we can conclude that without the information
from expression profiles, it becomes difficult for our model to learn unified
gene embeddings because of the decline in scores of different metrics.

Table 8: Comparison between with feature design and without feature design. Methods | ASW | AUC | iLISI | GC | CGR | NO
---|---|---|---|---|---|---
with feature | 0.75$\pm$0.01 | 0.78$\pm$0.02 | 0.53$\pm$0.01 | 0.73$\pm$0.04 | 0.65$\pm$0.02 | 0.31$\pm$0.00
w/o feature | 0.75$\pm$0.05 | 0.64$\pm$0.02 | 0.24$\pm$0.04 | 0.43$\pm$0.03 | 0.02$\pm$0.01 | 0.10$\pm$0.03
Figure 8: Gene embeddings based on different input information. (a) represents
the UMAPs of gene embeddings without expression profiles and (b) represents
the UMAPs of gene emebddings with expression profiles.

To demonstrate the necessity of Weighted Similarity Learning (WSL) design, we
design an experiment to analyze the performance of different model structures
with WSL or w/o WSL. According to Table 9, we found that MuSe-GNN with WSL
performed better than MuSe-GNN without WSL across all of the metrics in heart
tissue. Therefore, we need to utilize WSL manchism to effectively learn the
gene-gene relationships across different datasets.

Table 9: Comparison between with weighted similarity learning design and without weighted similarity learning design. Methods | ASW | AUC | iLISI | GC | CGR | NO
---|---|---|---|---|---|---
with WSL | 0.75$\pm$0.01 | 0.78$\pm$0.02 | 0.53$\pm$0.01 | 0.73$\pm$0.04 | 0.65$\pm$0.02 | 0.31$\pm$0.00
w/o WSL | 0.69$\pm$0.02 | 0.66$\pm$0.03 | 0.46$\pm$0.03 | 0.57$\pm$0.07 | 0.48$\pm$0.05 | 0.29$\pm$0.00

####  E.2 Hyper-parameter tuning

Hyper-parameter settings for various models can be divided into two groups: 1.
Methods such as PCA, Gene2vec, GIANT, WSMAE, GAE, VGAE, and MAE, which are
free of pre-training, are configured with the best hyper-parameters based on
parameter searching to ensure fairness. Details are included in Table 10. 2.
Method such as scBERT uses the pre-trained model. The optimal hyper-parameters
for MuSe-GNN are provided in Table 11. Number of trainable parameters for
different models is summarized in Table 12.

Table 10: Hyper-parameter candidates for benchmarking methods. Models | Hyper-parameters
---|---
PCA | Dim:[32,128]
Gene2vec | Dim:[32,128]
GIANT | Dim:[32,128]
WSMAE | LR:[1e-4,1e-2]; Batch size:[1000,3000]; Dim:[32,128]
GAE | LR:[1e-4,1e-2]; Batch size:[1000,3000]; Dim:[32,128]
VGAE | LR:[1e-4,1e-2]; Batch size:[1000,3000]; Dim:[32,128]
MAE | LR:[1e-4,1e-2]; Batch size:[1000,3000]; Dim:[32,128]
Table 11: Hyper-parameters list for MuSe-GNN. Parameter | Value |  |  |
---|---|---|---|---
Epoch | 2000 |  |  |
LR of encoder | 1e-4 |  |  |
LR of decoder | 1e-3 |  |  |
$\lambda_{c}$ | 1e-2 |  |  |
Dim of embeddings | 32 |  |  |
Sample size | 100 |  |  |
Table 12: Number of trainable parameters for different models Models | # of parameters (M)
---|---
MuSe-GNN | 349
GIANT | 260
VGAE | 210
WSMAE | 105
GAE | 52.5
MAE | 52.5

The average rank results for various hyper-parameters are presented in Table
13. We set the initial parameter settings as follows: Epoch = 2000, LR of
encoder = 1e-2, LR of decoder = 1e-3, $\lambda_{c}$ = 1e-2, Dim = 32, and
Sample size (for contrastive learning) = 100. Notably, for the value choice of
$\lambda_{c}$, we have two candidates with the same average rank result. To
further compare the gene embeddings generated by models based on these two
choices, we plot the UMAPs for these two models. According to Figure 9, we
observe that choosing $\lambda_{c}=0.1$ fails to integrate genes from dataset
G19. Therefore, we ultimately set $\lambda_{c}$ = 1e-2.

Additionally, the dimensions of gene embeddings should not be set to a very
large value, not only due to the results in our table but also because of
theoretical constraints and practical feasibility. Generally, an excessively
large number of embedding settings may cause the number of hidden layer
neurons to exceed the number of features in the input partial dataset,
increasing the probability of overfitting by introducing more noise [57].
Furthermore, a larger number of dimensions will consume more computational
resources, potentially leading to out-of-memory issues, such as when setting
the Dim of embeddings = 256.

Table 13: Hyper-parameter searching results. LR for encoder | Avg Rank | InfoNCE penalty | Avg Rank
---|---|---|---
0.01 | 3.17 | 1 | 3.5
0.001 | 2.50 | 0.1 | 2.13
0.0005 | 2.33 | 0.01 | 2.13
0.0001 | 1.83 | 0.001 | 2.25
LR for decoder | Avg Rank | Dim | Avg Rank
0.01 | 3.43 | 32 | 1.86
0.001 | 1.86 | 64 | 2.00
0.0005 | 2.29 | 128 | 2.14
0.0001 | 2.43 | 256 | 4.00
batch size | Avg Rank | sample size | Avg Rank
1000 | 3 | 50 | 2.38
1500 | 3 | 100 | 1.75
2000 | 1.5 | 150 | 2.88
3000 | 2.5 | 200 | 2.5
Figure 9: Gene embeddings on different $\lambda_{c}$.

####  E.3 Details about benchmarking experiments

The average ranks result for all the scRNA-seq datasets is displayed in Table
14. The experimental results for all the scRNA-seq datasets are displayed in
Tables 15 to 24 (unscaled). Based on the benchmark results for different
tissues, we can conclude that MuSe-GNN consistently outperforms other methods
according to specific metrics, such as iLISI and NO. Furthermore, the gene
embeddings learned by MuSe-GNN are relatively more stable than those from
other methods, as evidenced by the analysis of standard deviation. Considering
the overall performance, MuSe-GNN, as the top-performing method, also
surpasses the second-best methods from 20.1% to 97.5% and surpasses the third-
best methods from 26.9% to 99.7%.

Table 14: The average ranks for gene embeddings benchmarking across different tissues. Methods | Heart | Lung | Liver | Kidney | Thymus | Spleen | Pancreas | Cerebrum | Cerebellum | PBMC
---|---|---|---|---|---|---|---|---|---|---
PCA | 4.33 | 4.33 | 4.67 | 5.50 | 4.00 | 4.17 | 4.83 | 4.50 | 5.50 | 4.00
Gene2vec | 5.67 | 6.67 | 6.50 | 6.67 | 7.33 | 7.17 | 6.83 | 6.67 | 6.00 | 7.33
GIANT | 6.00 | 4.67 | 5.33 | 4.67 | 3.67 | 5.67 | 3.83 | 5.50 | 6.00 | 6.17
WSMAE | 5.83 | 5.33 | 4.83 | 4.83 | 3.83 | 5.17 | 5.00 | 5.33 | 4.33 | 5.00
GAE | 4.33 | 5.83 | 5.33 | 5.50 | 5.17 | 5.33 | 4.50 | 6.00 | 5.00 | 5.33
VGAE | 3.50 | 6.33 | 6.33 | 6.67 | 4.83 | 6.50 | 6.17 | 6.67 | 6.83 | 5.83
MAE | 7.00 | 5.17 | 5.17 | 5.17 | 5.50 | 5.17 | 5.67 | 5.33 | 5.33 | 5.33
scBERT | 5.67 | 4.50 | 4.33 | 3.17 | 7.83 | 3.83 | 4.83 | 3.17 | 3.17 | 3.67
MuSeGNN | 2.67 | 2.17 | 2.50 | 2.83 | 2.67 | 2.00 | 2.83 | 1.67 | 2.67 | 2.33
Table 15: Benchmark score table for Heart. Methods | ASW | AUC | iLISI | GC | CGR | NO
---|---|---|---|---|---|---
PCA | 0.77$\pm$0.01 | 0.64$\pm$0.00 | 0.43$\pm$0.00 | 0.49$\pm$0.02 | 0.30$\pm$0.01 | 0.27$\pm$0.00
Gene2vec | 0.82$\pm$0.03 | 0.76$\pm$0.00 | 0.09$\pm$0.01 | 0.65$\pm$0.03 | 0.00$\pm$0.00 | 0.10$\pm$0.02
GIANT | 0.84$\pm$0.01 | 0.55$\pm$0.00 | 0.26$\pm$0.02 | 0.48$\pm$0.03 | 0.10$\pm$0.01 | 0.01$\pm$0.00
WSMAE | 0.79$\pm$0.01 | 0.72$\pm$0.00 | 0.34$\pm$0.01 | 0.43$\pm$0.04 | 0.11$\pm$0.02 | 0.26$\pm$0.01
GAE | 0.79$\pm$0.01 | 0.98$\pm$0.00 | 0.38$\pm$0.01 | 0.45$\pm$0.02 | 0.12$\pm$0.02 | 0.28$\pm$0.01
VGAE | 0.80$\pm$0.01 | 0.97$\pm$0.00 | 0.39$\pm$0.01 | 0.44$\pm$0.02 | 0.18$\pm$0.03 | 0.28$\pm$0.00
MAE | 0.82$\pm$0.01 | 0.52$\pm$0.00 | 0.30$\pm$0.02 | 0.24$\pm$0.05 | 0.03$\pm$0.00 | 0.22$\pm$0.01
scBERT | 0.77$\pm$0.01 | 0.47$\pm$0.00 | 0.39$\pm$0.01 | 0.37$\pm$0.03 | 0.26$\pm$0.01 | 0.28$\pm$0.00
MuSe-GNN | 0.75$\pm$0.01 | 0.78$\pm$0.02 | 0.53$\pm$0.01 | 0.73$\pm$0.04 | 0.65$\pm$0.01 | 0.31$\pm$0.00
Table 16: Benchmark score table for Lung. Methods | ASW | AUC | iLISI | GC | CGR | NO
---|---|---|---|---|---|---
PCA | 0.79$\pm$0.00 | 0.69$\pm$0.00 | 0.58$\pm$0.00 | 0.70$\pm$0.01 | 0.07$\pm$0.01 | 0.10$\pm$0.00
Gene2vec | 0.71$\pm$0.40 | 0.76$\pm$0.00 | 0.00$\pm$0.00 | 0.89$\pm$0.01 | 0.00$\pm$0.00 | 0.01$\pm$0.00
GIANT | 0.90$\pm$0.01 | 0.49$\pm$0.00 | 0.12$\pm$0.03 | 0.88$\pm$0.03 | 0.07$\pm$0.01 | 0.01$\pm$0.00
WSMAE | 0.76$\pm$0.01 | 0.69$\pm$0.01 | 0.45$\pm$0.03 | 0.77$\pm$0.03 | 0.05$\pm$0.01 | 0.09$\pm$0.01
GAE | 0.73$\pm$0.03 | 0.86$\pm$0.00 | 0.23$\pm$0.04 | 0.78$\pm$0.03 | 0.02$\pm$0.00 | 0.07$\pm$0.01
VGAE | 0.51$\pm$0.10 | 0.87$\pm$0.01 | 0.00$\pm$0.00 | 0.87$\pm$0.03 | 0.00$\pm$0.00 | 0.01$\pm$0.01
MAE | 0.85$\pm$0.02 | 0.62$\pm$0.01 | 0.35$\pm$0.05 | 0.81$\pm$0.03 | 0.02$\pm$0.01 | 0.08$\pm$0.01
scBERT | 0.87$\pm$0.01 | 0.48$\pm$0.00 | 0.29$\pm$0.02 | 0.86$\pm$0.05 | 0.06$\pm$0.01 | 0.10$\pm$0.00
MuSe-GNN | 0.84$\pm$0.02 | 0.86$\pm$0.01 | 0.64$\pm$0.04 | 0.87$\pm$0.03 | 0.37$\pm$0.03 | 0.19$\pm$0.00
Table 17: Benchmark score table for Liver. Methods | ASW | AUC | iLISI | GC | CGR | NO
---|---|---|---|---|---|---
PCA | 0.77$\pm$0.00 | 0.58$\pm$0.00 | 0.39$\pm$0.00 | 0.69$\pm$0.01 | 0.21$\pm$0.00 | 0.15$\pm$0.00
Gene2vec | 0.63$\pm$0.34 | 0.69$\pm$0.00 | 0.00$\pm$0.00 | 0.88$\pm$0.01 | 0.00$\pm$0.00 | 0.03$\pm$0.01
GIANT | 0.90$\pm$0.01 | 0.55$\pm$0.00 | 0.40$\pm$0.02 | 0.44$\pm$0.03 | 0.11$\pm$0.01 | 0.00$\pm$0.00
WSMAE | 0.80$\pm$0.01 | 0.62$\pm$0.01 | 0.26$\pm$0.03 | 0.83$\pm$0.03 | 0.09$\pm$0.01 | 0.13$\pm$0.01
GAE | 0.74$\pm$0.02 | 0.84$\pm$0.00 | 0.17$\pm$0.02 | 0.80$\pm$0.01 | 0.05$\pm$0.01 | 0.14$\pm$0.01
VGAE | 0.49$\pm$0.05 | 0.87$\pm$0.00 | 0.00$\pm$0.00 | 0.87$\pm$0.02 | 0.01$\pm$0.00 | 0.01$\pm$0.00
MAE | 0.83$\pm$0.01 | 0.58$\pm$0.01 | 0.25$\pm$0.04 | 0.88$\pm$0.04 | 0.04$\pm$0.01 | 0.09$\pm$0.01
scBERT | 0.87$\pm$0.01 | 0.50$\pm$0.00 | 0.23$\pm$0.05 | 0.90$\pm$0.03 | 0.11$\pm$0.02 | 0.14$\pm$0.00
MuSe-GNN | 0.86$\pm$0.01 | 0.78$\pm$0.03 | 0.65$\pm$0.02 | 0.83$\pm$0.04 | 0.46$\pm$0.02 | 0.18$\pm$0.00
Table 18: Benchmark score table for Kidney. Methods | ASW | AUC | iLISI | GC | CGR | NO
---|---|---|---|---|---|---
PCA | 0.69$\pm$0.01 | 0.55$\pm$0.00 | 0.39$\pm$0.01 | 0.98$\pm$0.00 | 0.00$\pm$0.00 | 0.21$\pm$0.00
Gene2vec | 0.33$\pm$0.42 | 0.74$\pm$0.00 | 0.00$\pm$0.00 | 1.00$\pm$0.01 | 0.00$\pm$0.00 | 0.02$\pm$0.03
GIANT | 0.92$\pm$0.00 | 0.55$\pm$0.00 | 0.89$\pm$0.01 | 0.78$\pm$0.03 | 0.02$\pm$0.00 | 0.01$\pm$0.00
WSMAE | 0.74$\pm$0.05 | 0.68$\pm$0.01 | 0.08$\pm$0.04 | 0.99$\pm$0.00 | 0.01$\pm$0.00 | 0.15$\pm$0.04
GAE | 0.52$\pm$0.10 | 0.86$\pm$0.00 | 0.01$\pm$0.01 | 0.98$\pm$0.01 | 0.00$\pm$0.00 | 0.05$\pm$0.02
VGAE | 0.51$\pm$0.16 | 0.88$\pm$0.00 | 0.00$\pm$0.00 | 0.98$\pm$0.01 | 0.00$\pm$0.00 | 0.00$\pm$0.01
MAE | 0.84$\pm$0.04 | 0.52$\pm$0.01 | 0.17$\pm$0.07 | 1.00$\pm$0.00 | 0.00$\pm$0.00 | 0.16$\pm$0.04
scBERT | 0.89$\pm$0.02 | 0.47$\pm$0.00 | 0.59$\pm$0.18 | 1.00$\pm$0.00 | 0.05$\pm$0.01 | 0.30$\pm$0.00
MuSe-GNN | 0.88$\pm$0.03 | 0.71$\pm$0.02 | 0.84$\pm$0.04 | 0.98$\pm$0.02 | 0.31$\pm$0.04 | 0.31$\pm$0.01
Table 19: Benchmark score table for Thymus. Methods | ASW | AUC | iLISI | GC | CGR | NO
---|---|---|---|---|---|---
PCA | 0.76$\pm$0.01 | 0.79$\pm$0.00 | 0.93$\pm$0.00 | 0.44$\pm$0.05 | 0.00$\pm$0.00 | 0.00$\pm$0.00
Gene2vec | 0.00$\pm$0.00 | 0.75$\pm$0.00 | 0.00$\pm$0.00 | 0.58$\pm$0.03 | 0.00$\pm$0.00 | 0.00$\pm$0.00
GIANT | 0.88$\pm$0.02 | 0.55$\pm$0.00 | 0.98$\pm$0.01 | 0.80$\pm$0.02 | 0.00$\pm$0.00 | 0.00$\pm$0.00
WSMAE | 0.81$\pm$0.03 | 0.74$\pm$0.01 | 0.93$\pm$0.06 | 0.48$\pm$0.06 | 0.00$\pm$0.00 | 0.00$\pm$0.00
GAE | 0.72$\pm$0.12 | 0.81$\pm$0.00 | 0.89$\pm$0.10 | 0.69$\pm$0.05 | 0.00$\pm$0.00 | 0.00$\pm$0.00
VGAE | 0.77$\pm$0.05 | 0.84$\pm$0.01 | 0.89$\pm$0.05 | 0.44$\pm$0.10 | 0.00$\pm$0.00 | 0.00$\pm$0.00
MAE | 0.77$\pm$0.28 | 0.59$\pm$0.01 | 0.94$\pm$0.05 | 0.33$\pm$0.05 | 0.00$\pm$0.00 | 0.00$\pm$0.00
scBERT | 0.09$\pm$0.27 | 0.50$\pm$0.00 | 0.41$\pm$0.18 | 0.42$\pm$0.04 | 0.00$\pm$0.00 | 0.00$\pm$0.00
MuSe-GNN | 0.71$\pm$0.05 | 0.80$\pm$0.02 | 0.98$\pm$0.01 | 0.65$\pm$0.07 | 0.01$\pm$0.00 | 0.00$\pm$0.00
Table 20: Benchmark score table for Spleen. Methods | ASW | AUC | iLISI | GC | CGR | NO
---|---|---|---|---|---|---
PCA | 0.81$\pm$0.01 | 0.56$\pm$0.00 | 0.47$\pm$0.00 | 0.73$\pm$0.01 | 0.23$\pm$0.01 | 0.17$\pm$0.00
Gene2vec | 0.69$\pm$0.36 | 0.68$\pm$0.01 | 0.00$\pm$0.00 | 0.85$\pm$0.01 | 0.00$\pm$0.00 | 0.01$\pm$0.01
GIANT | 0.88$\pm$0.01 | 0.55$\pm$0.00 | 0.42$\pm$0.02 | 0.45$\pm$0.02 | 0.10$\pm$0.00 | 0.01$\pm$0.00
WSMAE | 0.79$\pm$0.01 | 0.61$\pm$0.01 | 0.26$\pm$0.03 | 0.86$\pm$0.02 | 0.08$\pm$0.01 | 0.15$\pm$0.01
GAE | 0.75$\pm$0.02 | 0.84$\pm$0.00 | 0.18$\pm$0.02 | 0.81$\pm$0.01 | 0.05$\pm$0.01 | 0.16$\pm$0.01
VGAE | 0.52$\pm$0.05 | 0.86$\pm$0.00 | 0.00$\pm$0.00 | 0.85$\pm$0.03 | 0.00$\pm$0.00 | 0.01$\pm$0.00
MAE | 0.84$\pm$0.02 | 0.56$\pm$0.01 | 0.30$\pm$0.04 | 0.92$\pm$0.02 | 0.03$\pm$0.00 | 0.12$\pm$0.00
scBERT | 0.87$\pm$0.01 | 0.50$\pm$0.00 | 0.34$\pm$0.06 | 0.94$\pm$0.02 | 0.11$\pm$0.01 | 0.16$\pm$0.00
MuSe-GNN | 0.86$\pm$0.01 | 0.79$\pm$0.02 | 0.70$\pm$0.02 | 0.89$\pm$0.03 | 0.47$\pm$0.01 | 0.19$\pm$0.00
Table 21: Benchmark score table for Pancreas. Methods | ASW | AUC | iLISI | GC | CGR | NO
---|---|---|---|---|---|---
PCA | 0.61$\pm$0.01 | 0.62$\pm$0.00 | 0.64$\pm$0.00 | 0.57$\pm$0.01 | 0.00$\pm$0.00 | 0.01$\pm$0.00
Gene2vec | 0.17$\pm$0.36 | 0.73$\pm$0.00 | 0.00$\pm$0.00 | 0.51$\pm$0.05 | 0.00$\pm$0.00 | 0.00$\pm$0.00
GIANT | 0.91$\pm$0.01 | 0.55$\pm$0.00 | 0.97$\pm$0.03 | 0.65$\pm$0.05 | 0.00$\pm$0.00 | 0.00$\pm$0.00
WSMAE | 0.67$\pm$0.09 | 0.70$\pm$0.01 | 0.66$\pm$0.10 | 0.49$\pm$0.06 | 0.00$\pm$0.00 | 0.00$\pm$0.00
GAE | 0.09$\pm$0.27 | 0.82$\pm$0.00 | 0.90$\pm$0.07 | 0.65$\pm$0.03 | 0.00$\pm$0.00 | 0.00$\pm$0.00
VGAE | 0.37$\pm$0.40 | 0.83$\pm$0.01 | 0.13$\pm$0.10 | 0.51$\pm$0.04 | 0.00$\pm$0.00 | 0.00$\pm$0.00
MAE | 0.71$\pm$0.11 | 0.56$\pm$0.01 | 0.83$\pm$0.07 | 0.37$\pm$0.05 | 0.00$\pm$0.00 | 0.00$\pm$0.00
scBERT | 0.76$\pm$0.05 | 0.47$\pm$0.00 | 0.54$\pm$0.15 | 0.39$\pm$0.04 | 0.00$\pm$0.00 | 0.01$\pm$0.00
MuSe-GNN | 0.72$\pm$0.08 | 0.65$\pm$0.01 | 0.95$\pm$0.13 | 0.62$\pm$0.06 | 0.01$\pm$0.00 | 0.00$\pm$0.00
Table 22: Benchmark score table for Cerebrum. Methods | ASW | AUC | iLISI | GC | CGR | NO
---|---|---|---|---|---|---
PCA | 0.82$\pm$0.01 | 0.63$\pm$0.00 | 0.34$\pm$0.00 | 0.95$\pm$0.00 | 0.12$\pm$0.00 | 0.18$\pm$0.00
Gene2vec | 0.07$\pm$0.23 | 0.69$\pm$0.00 | 0.00$\pm$0.00 | 0.99$\pm$0.00 | 0.00$\pm$0.00 | 0.01$\pm$0.01
GIANT | 0.89$\pm$0.00 | 0.55$\pm$0.00 | 0.40$\pm$0.02 | 0.51$\pm$0.02 | 0.03$\pm$0.00 | 0.01$\pm$0.00
WSMAE | 0.80$\pm$0.02 | 0.67$\pm$0.01 | 0.21$\pm$0.04 | 0.97$\pm$0.01 | 0.03$\pm$0.01 | 0.13$\pm$0.02
GAE | 0.68$\pm$0.04 | 0.83$\pm$0.00 | 0.03$\pm$0.01 | 0.95$\pm$0.01 | 0.02$\pm$0.00 | 0.13$\pm$0.01
VGAE | 0.47$\pm$0.09 | 0.84$\pm$0.01 | 0.00$\pm$0.00 | 0.98$\pm$0.01 | 0.00$\pm$0.00 | 0.01$\pm$0.00
MAE | 0.83$\pm$0.02 | 0.53$\pm$0.01 | 0.20$\pm$0.05 | 0.99$\pm$0.01 | 0.02$\pm$0.01 | 0.14$\pm$0.01
scBERT | 0.89$\pm$0.01 | 0.48$\pm$0.00 | 0.50$\pm$0.07 | 0.99$\pm$0.00 | 0.10$\pm$0.01 | 0.17$\pm$0.00
MuSe-GNN | 0.90$\pm$0.01 | 0.73$\pm$0.02 | 0.79$\pm$0.03 | 0.99$\pm$0.01 | 0.54$\pm$0.01 | 0.21$\pm$0.00
Table 23: Benchmark score table for Cerebellum. Methods | ASW | AUC | iLISI | GC | CGR | NO
---|---|---|---|---|---|---
PCA | 0.82$\pm$0.00 | 0.63$\pm$0.00 | 0.34$\pm$0.00 | 0.95$\pm$0.00 | 0.12$\pm$0.01 | 0.18$\pm$0.00
Gene2vec | 0.31$\pm$0.40 | 0.71$\pm$0.00 | 0.00$\pm$0.00 | 1.00$\pm$0.00 | 0.00$\pm$0.00 | 0.02$\pm$0.02
GIANT | 0.89$\pm$0.00 | 0.55$\pm$0.00 | 0.34$\pm$0.03 | 0.73$\pm$0.02 | 0.03$\pm$0.00 | 0.01$\pm$0.00
WSMAE | 0.83$\pm$0.02 | 0.66$\pm$0.01 | 0.26$\pm$0.04 | 1.00$\pm$0.00 | 0.06$\pm$0.01 | 0.29$\pm$0.02
GAE | 0.74$\pm$0.02 | 0.83$\pm$0.00 | 0.11$\pm$0.08 | 0.99$\pm$0.01 | 0.07$\pm$0.04 | 0.25$\pm$0.04
VGAE | 0.19$\pm$0.26 | 0.83$\pm$0.01 | 0.00$\pm$0.00 | 0.99$\pm$0.01 | 0.00$\pm$0.00 | 0.00$\pm$0.00
MAE | 0.86$\pm$0.01 | 0.51$\pm$0.00 | 0.16$\pm$0.03 | 1.00$\pm$0.00 | 0.02$\pm$0.00 | 0.26$\pm$0.03
scBERT | 0.90$\pm$0.01 | 0.47$\pm$0.00 | 0.54$\pm$0.06 | 1.00$\pm$0.00 | 0.15$\pm$0.02 | 0.31$\pm$0.01
MuSe-GNN | 0.92$\pm$0.01 | 0.64$\pm$0.02 | 0.86$\pm$0.02 | 0.98$\pm$0.02 | 0.65$\pm$0.00 | 0.38$\pm$0.00
Table 24: Benchmark score table for PBMC. Methods | ASW | AUC | iLISI | GC | CGR | NO
---|---|---|---|---|---|---
PCA | 0.86$\pm$0.01 | 0.53$\pm$0.00 | 0.46$\pm$0.00 | 0.66$\pm$0.01 | 0.22$\pm$0.01 | 0.13$\pm$0.00
GENE2VEC | 0.00$\pm$0.00 | 0.59$\pm$0.00 | 0.00$\pm$0.00 | 0.92$\pm$0.01 | 0.00$\pm$0.00 | 0.00$\pm$0.00
GIANT | 0.85$\pm$0.01 | 0.55$\pm$0.00 | 0.38$\pm$0.02 | 0.55$\pm$0.03 | 0.04$\pm$0.00 | 0.00$\pm$0.00
WSMAE | 0.81$\pm$0.01 | 0.58$\pm$0.01 | 0.29$\pm$0.02 | 0.81$\pm$0.05 | 0.10$\pm$0.02 | 0.12$\pm$0.01
GAE | 0.74$\pm$0.01 | 0.82$\pm$0.00 | 0.15$\pm$0.03 | 0.79$\pm$0.03 | 0.06$\pm$0.02 | 0.12$\pm$0.01
VGAE | 0.67$\pm$0.04 | 0.89$\pm$0.00 | 0.00$\pm$0.00 | 0.93$\pm$0.02 | 0.02$\pm$0.00 | 0.01$\pm$0.00
MAE | 0.85$\pm$0.02 | 0.53$\pm$0.01 | 0.30$\pm$0.05 | 0.94$\pm$0.03 | 0.02$\pm$0.00 | 0.10$\pm$0.01
scBERT | 0.88$\pm$0.01 | 0.52$\pm$0.00 | 0.48$\pm$0.05 | 0.93$\pm$0.03 | 0.15$\pm$0.02 | 0.11$\pm$0.00
MuSe-GNN | 0.86$\pm$0.01 | 0.79$\pm$0.02 | 0.75$\pm$0.02 | 0.88$\pm$0.02 | 0.54$\pm$0.01 | 0.18$\pm$0.00

####  E.4 Details about function cluster identification

In this paper, we utilized Leiden cluster method [91] to identify genes with
common functions based on the value of their embeddings. The genes with
distance smaller than a community-discovering based threshold (known as co-
embedded genes) will be identified into a group. Here we sample 100 genes in
our final gene embedding result from heart tissue to illustrate the
distribution of such distance.

Figure 10: Heatmap of gene embedding distances. The heatmap color represents
the distance between genes, with brighter colors indicating greater distances
and more distinct functions.

Based on Figure 10, we can conclude that our method has the ability to
identify genes with similar embeddings, and further, it can be used to derive
genes with similar functions.

All of the experiments of this paper are finished based on Intel Xeon Gold
6240 Processor and one NVIDIA A100 GPU. The RAM setting is 30 GB.

###  Appendix F Extra Discussion about Model Design

This section aims to demonstrate the reasoning behind our model design in
greater details. As a novel graph neural network model, MuSe-GNN offers three
significant advantages.

Firstly, our method’s graph construction process is robust and reliable. Graph
construction is a critical step for GNN models [117], and we place
considerable emphasis on both the graph construction method and data quality
throughout our work. In single-cell data analysis research, some researchers
model cell-cell similarity based on the nearest neighbor graph [99, 94]. This
approach may be affected by batch effects resulting from different sequencing
profiles, leading to incorrect conclusions. For instance, in scRNA-seq data
imputation tasks, such a design may result in a high false positive rate [5].
A gene co-expression network derived from a robust statistical framework is
more trustworthy, and biological experiments can validate the co-expression
performance of various genes. Moreover, to understand the difference of co-
expression networks from different tissues or omics data, we calculate the
Edge Signal-Noise Ratio (ESNR) [22] for each network, and figure out that the
co-expression network based on scRNA-seq has the best quality, while the co-
expression networks from other omics have very small ESNR. Such conclusions
are also consistent with the model performance.

Secondly, our model utilizes a sophisticated and stable training strategy. By
employing the cosine similarity penalty, we ensure that genes with the same
functions from different tissues or techniques can be co-embedded in similar
positions. Additionally, we use graph auto-encoder and contrastive learning
strategies to retain the original biological information of various datasets.
Each of our training strategies serves a specific purpose.

Finally, our method showcases the high effectiveness of MMML in addressing
heterogeneous data learning challenges. MMML can successfully learn similar
information across multimodal data and represent it in a unified space through
model fusion [55]. This approach can be applied to more complex data studies
in the future, such as integrating data from Magnetic Resonance Imaging (MRI)
results, spatial resolution image information, and current expression profiles
into a single space. Furthermore, GNN remains the primary model for handling
tasks related to graph-structured data. By combining these methods and
uncovering biological insights, our research contributes to a deeper
understanding of efficient human genome representation.

###  Appendix G Selecting anchors of common functional genes

In this paper, we employ explicit common functional genes as anchors across
different datasets to implement regularization in our training process. The
common functional genes we select can potentially recover a portion of the
functional overlap of the genes, while the remaining implicit genes can be
recovered during the training process.

In this paper, we choose to use common highly variable genes (HVGs) combined
from different datasets as anchors, and employ their co-expression network
overlap as weights. An example is illustrated in Figure 12 case 1. In addition
to the proof provided by the ablation test results, there are four reasons for
this choice: 1. HVGs represent the active genes of a expression profile
(dataset), which are correlated to gene functions. 2. For datasets from the
same tissue, the overlap of HVGs and the overlap of neighbor genes (known as
the co-expressed genes in the same graph) of HVGs pair are similar. This
demonstrates a type of consistency in biological function inference. 3. For
datasets from different modalities (tissues, techniques, etc.), only using the
overlap of HVGs is not reliable. The overlap of neighbor genes of HVGs pair
from different tissues shows functional similarity based on hierarchical
clustering. We can observe the clear difference of HVGs overlap score and
neighbors overlap score of HVGs by comparing Figure 11 (a) to (b).
Furthermore, we can consider an extreme case, illustrated in Figure 12 case 2.
If two HVGs from different tissues do not contain shared neighbors, the
training process will be processed for such genes without weighted
constraints. 4. Since it is quite challenging to select all the true genes
with similar functions, we can perform a trade-off between learning the
similarity of potential functionally similar genes and learning the difference
of potential functionally distinct genes. The neighbor genes overlap in common
HVGs is significantly different from the overlap in different HVGs for graphs
from the same tissue, as shown in Figure 13. Such statistical results also
support our choice of using common HVGs as anchors.

Figure 11: Cluster heatmaps of anchor genes for multimodal data. (a)
represents the cluster heatmap based on HVGs overlap, (b) represents the
cluster heatmap based on the neighbor genes overlap of HVGs.  Figure 12: The
need for weight similarity learning illustrated by two cases.  Figure 13:
Boxplots for the neighbor overlap score of HVGs from different tissues. We
used Kolmogorov–Smirnov test [97] to compare the difference between the two
groups of genes. The significance level of group difference is shown by the
stars ($***:\text{P-value}\leq 0.01$). (a) represent the distribution of
overlap scores for heart tissue, grouped by common HVGs and different HVGs.
(b) represent the distribution of overlap score for cerebrum tissue, grouped
by common HVGs and different HVGs.

In conclusion, we aim to generate similar embeddings for functionally related
genes, incorporating common HVGs and their neighbors in co-expression
networks, by utilizing the regularization term associated with common HVGs.
Conversely, embeddings for functionally distinct genes should differ, enabling
their identification through community detection-based clustering algorithms.

###  Appendix H Metrics details

Here we assume that HVGs from the same tissues with the same name have similar
functions across different datasets, and we also consider evaluating the
preservation of the known gene co-expression relation by using node
similarity. These metrics include:

  * •

edge AUC: Edge Area under the ROC Curve (AUC) is a metric widely used in the
edge prediction task. We calculate the AUC score between the multiplication of
normalized embeddings with Sigmoid threshold ($Sigmoid(zz^{T})$) and the true
edge relation. This metric can reflect the co-expression information we have
in our embedding space.

  * •

common genes ASW: Common genes ASW represents the ASW score we calculated
based on common genes across different tissues. Since the overlap genes of
different datasets are not very large, based on the assumption that different
highly variable genes have similar functions in different datasets, we can
utilize common genes ASW as a metric to evaluate our integration performance.

  * •

common genes Graph connectivity: Graph connectivity is a metric to evaluate
the connectivity based on the ratio between the number of nodes in the largest
connected component of the graph constructed by the gene cluster and the
number of nodes belonging to the graph. We calculate the graph connectivity
score for each gene cluster and take the average.

  * •

common genes iLISI: The inverse Simpson’s index (iLISI) determines the number
of common genes that can be drawn from a neighbor list before one dataset (or
datasets from the same tissue) is observed twice. Therefore, this metric can
be used to evaluate the integration level of the common genes across different
datasets or tissues based on the choice of our index type.

  * •

common genes ratio: The common genes ratio is defined based on the weighted
average of one minus the ratio of unique genes from one cluster to the total
genes in the cluster across all the clusters. We utilize the Leiden method to
generate the clusters. This metric can be used to evaluate the level of
integration of the same functional genes by a given method.

  * •

neighbors overlap: The neighbor overlap value is defined based on the
neighbors of one gene in the co-expression network. For each cluster, we
calculate the similarity of genes’ neighbors among different datasets and
perform a weighted average of the similarity. A larger value means that the
method is able to integrate similar genes.

###  Appendix I Shared transcription factors and pathways analysis

Transcription factors (TFs) are crucial in transcriptional regulation [37]. By
examining the overlap of TFs across diverse datasets, we can assess the
similarities between distinct modal data in terms of transcriptional
regulation. We employed the TRRUST v2 database [37] to identify regulators for
each gene within the gene embedding space. Subsequently, for every cluster, we
extracted the transcription factors associated with the genes and computed the
transcription factor overlap for the various multimodal data linked to the
genes. We iteratively repeated this process for each cluster and documented
the number of shared transcription factors based on different combinations of
biological multimodal data. The results are illustrated in Figure 14.

Figure 14: Results of common TFs discovery.

Figure 14 reveals that tissues with similar functions exhibit greater TF
overlap, such as the heart and lung group, and the spleen and liver group.
These findings demonstrate that our generated embeddings effectively capture
the similarities between various data types. ComplexUpset [51] was utilized to
create this figure.

Figure 15: The proportion of tissue-specific TFs enrichment in the target
tissue (blue) and in other tissues (orange). Here the significant enrichment
in blue groups are marked by stars ($*:\text{P-value}\leq
0.1;**:\text{P-value}\leq 0.05;***:\text{P-value}\leq 0.01;$).

Additionally, we can emphasize the enrichment status of tissue-specific TFs,
as depicted in Figure 15. We utilized the TF information from [107] to compute
the enrichment proportion for tissue-specific TFs enriched in the given tissue
and other tissues. A one-tailed Fisher’s exact test was also conducted to
assess the significance of enrichment differences. The figure suggests that
our gene embeddings consistently enrich tissue-specific TFs across all
evaluated tissues. Notably, for some tissues, the enrichment is significant,
which aligns with the conclusions from [16]. However, in other tissues, the
difference is not significant. A potential explanation is that such tissues
may play a critical role in the interaction between different tissues for
specific biological processes.

Figure 16: Top1 enriched pathways across different gene function clusters
grouped by major tissue in each cluster. Here the color of bubbles represents
the $-ln(FDR)$ value and the size of bubbles represents the value of rich
factor for different pathways.

In Figure 16, we determined the most enriched pathway for each cluster and
classified it into the most representative tissue within that cluster by GOEA.
The representative tissue was chosen based on the highest proportion of the
source tissue for genes. This figure allows us to conclude that the heart
tissue is the most versatile in the human body, followed by liver and lung
tissues, which are also crucial tissues in humans.

###  Appendix J Multi-species gene embeddings

Figure 17: Gene representation learning results for multimodal biological data
with different species. (a) represents the UMAPs of gene embeddings colored by
tissue type. (b) represents the UMAPs of gene embeddings colored by common
function groups. (c) is a Sankey plot [80] to show the genes overlap of
different species in the same clusters. (d) shows the top5 pathways related to
the genes in the special cluster discovered by GOEA. The bubble plots in this
paper were created based on ggplot2 [102].  Figure 18: Compare the proportion
of orthology genes discovered by MuSe-GNN and random selection.

Integrating datasets from different species based on cells is challenging
because we always need extra information, either orthology genes (genes from
different species whose evaluation is only related to speciation events)
information or protein embeddings information [76]. However, using genes as
anchors simplifies this task. We examined gene embeddings generated by MuSe-
GNN across three mammalian species, including humans, lemurs, and mice. As
shown in Figure 17 (a) and (b), our method effectively integrated information
from different species into a shared space, preserving gene similarity across
various techniques, tissues, and species. We also observed the co-embedding of
genes from different species but the same tissue, supporting the similarity of
evaluation direction. Since the enrichment of orthology genes can also be used
as proof to evaluate whether genes with similar functions are grouped [32], we
collected genes from different clusters with three species pairs and computed
the proportion of these genes in the orthology genes database [4].
Furthermore, we compared the results based on gene embeddings to random
selection as a null approach. From Figure 18, we could conclude that gene
embeddings from our model can effectively enrich genes from different species
but with similar functions.

Furthermore, Figure 17 (c) reveals that genes from different species are
distributed across multiple common function clusters. To explore shared gene
functions across species, we selected one common function cluster containing
genes from all three species for GOEA analysis. Since most genes in this
cluster originated from heart tissue, we anticipated pathway enrichment
related to multicellular systems and circulatory systems, as displayed in