GDAs-OT: A Prediction Method of Gene-Disease Associations Based on Optimal Transport for Identifying Genes Related to Immune-Related Adverse Events

1 Introduction

Cancer remains a critical threat to global public health characterized by its substantial incidence and mortality. In recent years, the paradigm of cancer treatment has been revolutionized by immunotherapy, with Immune Checkpoint Inhibitors (ICIs) emerging as a representative modality [Lee et al., 2016]. With inhibitory pathways, most notably the PD-1/PD-L1 and CTLA-4 pathways, ICIs can restore cytotoxic T-cell activity and enable effective antitumor immune responses [Postow et al., 2018]. However, this systemic immune activation is frequently accompanied by immune dysregulation. A growing body of clinical observations indicate that ICIs administration can provoke inflammatory toxicity in diverse organ systems, collectively termed immune-related Adverse Events (irAEs), which represents a major limitation to both treatment adherence and patient prognosis [RamosCasals et al., 2020]. Identifying genes associated with irAEs can help to further explore the underlying mechanisms of irAE occurrence, which in turn facilitates the refinement of safer immunotherapy protocols. In this work, we treat irAEs as diseases, and accordingly formulate the task as a GeneDisease Associations (GDAs) prediction problem, aiming to determine whether a given gene is associated with a specific disease.

∗Corresponding author

GDAs prediction is an important yet challenging task. Conventional experimental methods, such as linkage analysis [Ott et al., 2015] and Genome-Wide Association Studies (GWAS) [Manolio, 2010], are often expensive and require time-consuming "wet-lab" validation when applied to a large scale gene sets. To address these issues, a range of computational methods have been proposed based on conventional machine learning techniques [Natarajan and Dhillon, 2014; Isakov et al., 2017; Zeng et al., 2019] and graph neural networks [Li et al., 2023; Si et al., 2024; Meng et al., 2024], among which graph neural network-based approaches have shown particularly encouraging performance. Despite these advances, most existing graph neural network-based methods overlook disease-disease and disease-gene associations that have not yet been observed but biologically plausible. In this study, we refer to such missing links as false-negative associations. Concretely, the disease-disease associations used in training are incomplete. For gene-disease associations, confirmed links are treated as positive pairs, whereas links without established evidence are typically viewed as negative pairs. We examine the number of associations reported in different public releases of the gene-disease database DisGeNet [Pi˜nero et al., 2016; Pi˜nero et al., 2020; Pi˜nero et al., 2021] (see Table 1) and found that, many "negative" pairs in earlier versions are subsequently validated as positive associations in later releases. This temporal evolution of knowledge underscores the limitation of static negative sampling and highlights the need for uncertainty-aware or noiserobust graph inference. To address these issues, we propose GDAs-OT, a novel GDAs prediction method based on optimal transport. Optimal Transport (OT) was originally formulated as the problem of transforming one probability distribution into another at minimal cost, providing a principled measure of the overall structural discrepancy between distributions [Peyr´e et al., 2019].

We construct a heterogeneous network comprising a disease-disease graph, a gene-gene graph, and a gene-disease graph using publicly available data, with node features initialized using pre-trained PubMedBERT [Gu et al., 2021], a domain-specific transformer-based model pre-trained on biomedical corpora. In the disease-disease graph, our goal is to identify previously unrecorded but potentially meaningful positive pairs to enrich the disease association structure. Inspired by line graph theory [Harary and Norman, 1960], we compute edge features by combining the features of the two endpoint nodes and categorize edges into positive edges (with observed links) and unobserved edges (links that are currently absent). We then formulate an OT problem, minimizing the Wasserstein distance [Peyr´e et al., 2019] between the distributions of observed and unobserved edge features. Unobserved edges that receive transport mass exceeding a predefined threshold are identified as latent positives and used to enrich the graph's topological structure.

For the gene-disease graph, the focus shifts to selecting the most reliable negative pairs from unobserved links in order to reduce data contamination. Specifically, we minimize the Gromov-Wasserstein (GW) discrepancy [Peyr´e et al., 2016] to obtain an OT matrix between disease and gene nodes, and treat node pairs with the smallest transport mass as negative samples. After refining both the disease-disease and genedisease graphs, we map nodes from the disease-disease graph and the gene-gene graph into a shared embedding space using separate Graph Convolutional Network (GCN) encoders [Kipf and Welling, 2016]. Guided by the gene-disease relations, the learned embeddings pull positive gene-disease pairs closer while pushing negative pairs farther apart. Finally, the optimized node representations are fed into a multilayer perceptron (MLP) to produce the final GDA predictions.

We benchmark GDAs-OT against six state-of-the-art (SOTA) models for Gene-Disease Association tasks. Our framework achieves superior performance across all primary evaluation metrics, including Accuracy, AUC, F1-score, and AUPR. Ablation studies further confirm the contribution of each component in the framework. In addition, we curate several common irAEs and treat them as diseases, then use the trained GDAs-OT model to identify genes most strongly associated with each irAE. Existing evidence supports the reliability of the genes prioritized by GDAs-OT, and the model also highlights several potentially relevant genes that have not yet been reported, with meaningful implications for cancer immunotherapy. Our contributions can be concluded as:

• We proposed GDAs-OT, a novel prediction method of GDAs based on optimal transport for identifying genes related to immune-related adverse events, which was benchmarked against SOTA models, achieving the best performance.
• We partitioned edges in the disease-disease graph into observed positive edges and unobserved edges, learn an optimal transport plan by minimizing the Wasserstein distance, and use the resulting transport mass to augment the graph with high-confidence potential positive edges, for better capturing disease structural relationships.
• We learned an optimal transport matrix between disease nodes and gene nodes by minimizing the GromovWasserstein discrepancy, and leveraged it to construct a high-confidence set of negative pairs. These negatives guide representation learning in the embedding space and effectively mitigate performance degradation due to label noise and data contamination.
• We curateed several common irAEs as diseases and used GDAs-OT to prioritized associated genes, including plausible candidate genes that have not yet been reported, offering new directions for downstream biological investigation and cancer immunotherapy research.

2 RelatedWork

In the early stages, exploring the relationship between genes and diseases was more regarded as a priority ranking problem of disease genes. Specifically, given a set of known disease-associated genes and a set of candidate genes for a particular genetic disease, the objective of disease gene prioritization is to rank the candidate genes according to their potential relevance to the disease [Tran et al., 2020]. For example, [Natarajan and Dhillon, 2014] used the matrix completion method to integrate the features of multiple datasets to predict GDAs. [Isakov et al., 2017] proposed a machine learning-based gene priority ranking algorithm to infer risk genes for Inflammatory Bowel Disease (IBD). And [Tran et al., 2020] constructed a heterogeneous graph to fit a regularized linear SVM, thereby solving the gene priority ranking problem. XGDAG [Mastropietro et al., 2023] utilized a Markov process with restart to assign pseudo labels to unknown disease genes, and then used the filtered data to train GraphSAGE to achieve the ranking of candidate genes.

Recently, with the development of graph neural networks and the generation of a large amount of disease and gene data, exploring the relationship between genes and diseases has gradually been regarded as a link prediction problem. Specifically, it involves inferring the relationship between disease nodes and gene nodes. HOGCN [Kishan et al., 2021] implemented a GCN method that integrates high-order neighborhood information. Compared to only considering neighbor information, HOGCN achieved better results in GDAs prediction. GlaHGCL [Si et al., 2024] integrates miRNA, disease and gene data, combines global and local contrastive learning, and optimizes node embeddings in heterogeneous graphs for GDAs prediction. FusionGDA [Meng et al., 2024] proposed to use a pre-trained model to fuse protein sequences and disease descriptions to improve prediction performance.

3 Methodology

3.1 Preliminary

Biomedical Pre-Trained Language Models

These models are pre-trained on extensive biomedical datasets, enabling them to learn rich semantic representations from gene and disease texts. Once pre-training is completed, the models can be directly applied to obtain text embeddings. We compared BioGPT [Luo et al., 2022], LinkBERT [Yasunaga et al., 2022], and PubMedBERT in our experiments, and PubMedBERT was selected due to its superior performance. See the Supplementary Material Section 1 for details.

Graph Convolutional Network

GCN is a learning model for graph-structured data that learns node representations by aggregating information from neighboring nodes. Given a graph G(V,E), where V , E denote the sets of nodes and edges, respectively, let A[Aij]n×n be the adjacency matrix. If there exists an edge between node i, j, then Aij = 1; otherwise, Aij = 0. For the l - th layer, the node feature representation Hl is updated by aggregating features from the node itself and its first-order neighbors, yielding the (l +1)-th layer representation Hl+1 as follow:

where ˜ A = A + I denotes the adjacency matrix with added self-loops, ˜D is the corresponding degree matrix, W l represents the trainable parameters of the l - th layer, and σ(·) is a nonlinear activation function.

3.2 Data Preprocessing

To construct the gene-gene graph, disease-disease graph, and gene-disease graph required by GDAs-OT, we integrated data resources from multiple authoritative biomedical databases. Specifically, the associations between diseases and their textual descriptions originated from the Disease Ontology (DO) database [Schriml et al., 2022]. The relationships between genes originated from the HumanNet database [Kim et al., 2022], the functional description information of genes came from the Gene Ontology (GO) database [Ashburner et al., 2000], and the GDAs data were taken from the DisGeNet database [Pi˜nero et al., 2021]. Considering the differences in naming conventions and data coverage among different databases, certain filtering was performed. Based on disease names and UMLS Concept Unique Identifiers (CUIs), diseases in DisGeNet that were not present in DO were excluded. According to gene names, the gene associations that did not exist in HumanNet or GO in DisGeNet were also excluded. The processed data constitute the final dataset used in this paper's experiments. See Table 2.

After completing data preprocessing, we construct genegene graph Gg = G(V g ,Eg ) and disease-disease graph Gd = G(V d ,Ed ). Here, V g = {vg 1 , . . . , vg η } and V d = {vd 1 , . . . , vd ξ } denote the sets of gene nodes and disease nodes, respectively, while Eg and Ed represent the corresponding edge sets. Node features are initialized from textual descriptions of genes and diseases, which are encoded using the medical pre-trained model PubMedBERT. For any node v ∈ V , its feature representation is denoted as xv ∈ Rl . The feature matrix of all nodes is given by X ∈ R|V |×l .

3.3 Disease-Disease Graph Structure Expansion via Optimal Transport

In the disease graph Gd , to alleviate the false negative edge issue, we formulate the problem as an OT task between positive edges and unobserved edges defined on the same metric measure space. Specifically, positive edges are treated as the source distribution, while unobserved edges serve as the target distribution. By solving the minimization problem associated with the Wasserstein distance, we obtain an edge-level OT matrix, which is then used to structure the graph structure with potential positive associations. The set of positive edges consists of known disease-disease relationships. For node pairs that are not directly connected in the graph, we regard them as unobserved edges. Since the number of unobserved edges can be extremely large, we limit their scale by keeping only the five closest neighbors for each node according to cosine similarity of node features. Following the idea of line graph, we convert edges into nodes and derive edge representations accordingly. Let U = {ui|ei ∈ E} denote the set of edge features. For an edge ei = (vdw , vd v ∈ E), its feature is obtained by concatenating the features of its two endpoint nodes, which can be written as ui = [xdw , xdv ] ∈ R2d . We then divide U into a positive edge feature set U s ⊂ U and an unobserved edge feature set U t ⊂ U. Let the edge cost ma-

trix be denoted as C = [cij]m×n, where cij = 1- ui⊤·uj ||ui ||2||uj ||2 represents the cost of transporting unit mass from the i - th positive edge to the j - th unobserved edge. Here, m and n denote the numbers of edges in U s and U t , respectively. In the absence of prior knowledge, we assume a uniform mass distribution over edges, leading to the distributions on positive edges and unobserved edges as ps = 1m 1m, pt = 1n 1n. The Wasserstein distance between U s and U t is formulated as

Algorithm 1 Expansion of disease-disease graph structure

Input: G(V d ,Ed ), U s , U t , m = |U s |, threshold hyperparameter β.

Output: ˜G(V d , ˜Ed ).

• 1: Calculate the cost matrix C.
• 2: Initialize ˜Ed ←Ed .
• 3: Obtain T ∗ by solving Eq.(3).
• 4: Calculate τw ←min(1, 1 βm ).
• 5: for each unlabeled edge ej ∈ U t do
• 6: if maxi T ∗ ij > τw then
• 7: ˜Ed ← ˜Ed ∪ {ej}.
• 8: end if
• 9: end for
• 10: return ˜Gd (V d , ˜Ed ).

graph. These negative pairs will guide the representation learning process. Since this is an OT problem for two different metric measure spaces, we consider solving the OT matrix between genes and diseases by minimizing the GW discrepancy. The GW discrepancy can be viewed as a natural extension of the GW distance and is more suitable for machine learning applications [Peyr´e et al., 2016]. The dissimilarity between each internal node within a distribution is used to characterize the respective graph structure, thereby achieving the optimal structural alignment between two different metric measure spaces without relying on the comparability of crossdomain features. For Gg and ˜Gd , the intra-graph dissimilarity matrices are denoted as C g = [cgij ]η×η and C d = [cdi′ j ′ ]ξ×ξ,

mass of both gene and disease nodes are assumed to follow uniform distributions, i.e., pg = 1 η 1η, pd = 1ξ 1ξ. Following the common practice [M´emoli, 2011], we define the metric

where T = [Tij]m×n is the transport matrix. Each Tij corresponds to the amount of mass transported from a positive edge ei ∈ U s to an unobserved edge ej ∈ U t . < ·, · > denotes the inner product between matrices. ∏(ps, pt) = {T ∈ Rm×n |T 1n = ps, T ⊤1m = pt}.

This problem can be solved using the Sinkhorn-Knopp algorithm [Cuturi, 2013]. We consider the entropy regularization form of Eq.(2):

OT matrix T ∗ is calculated as

where diag(·) converts a vector to a diagonal matrix, K = exp (- Cλ ), λ = 0.1 is the regularization coefficient. Entropy regularization term H(T) = - ∑i,j Tij logTij. a ∈ Rm, b ∈ Rn is calculated through Sinkhorn iterations:

Thus, potential positive pairs in the unobserved edges is defined as an edge that receives the maximum transported mass from positive edges in the OT solution and whose value exceeds the threshold τw =min(1, 1 βm ), where β is a threshold weight. Specifically, for any ej ∈ U t , if maxi T ∗ ij > τw, it is treated as a positive pair. The disease-disease graph after edge completion is denoted as ˜Gd = ˜G(V d , ˜Ed ), with ˜Ed being the expanded set of positive disease-disease edges. The pseudo-code is given in Algorithm 1.

3.4 Negative Pairs Selection for the Gene-Disease Graph via Optimal Transport

Before conducting node representation learning, we need to select reliable negative pairs from the existing gene-disease measure spaces of Gg and ˜Gd as (C g , pg) ∈ Rη×η × ∑η and (C d , pd) ∈ Rξ×ξ × ∑ξ . The GW discrepancy between (C g , pg) and (C d , pd) is defined as

where ∏ (pg, pd) = ¯ T ∈ Rη×ξ | ¯ T1ξ = pg, ¯ T ⊤1η = pd}. ¯ T denotes the gene-disease transport matrix, and ¯ Tij represents the transported mass from gene node vi ∈ V g to disease node vj ∈ V d . L(C g ,C d , ¯ T) = [Ljj ′ ] ∈ R|V g |×|V d | and each Ljj ′ = ∑i,i′ L(cgij , cdi′ j ′) ¯ Tii′ , L(·, ·). The GW discrepancy allows flexible choices of loss functions, such as the squared loss L(a, b) = 12 |a - b|2 or the Kullback-Leibler divergence L(a, b) = a log ab - a + b. After experimental analysis(see Supplementary Material Section 2), we adopt the squared loss in this work.

To efficiently solve Eq.(6), we introduce an entropic regularization term and reformulate it as

Following the method proposed by [Peyr´e et al., 2016], this non-convex optimization problem is solved using projected gradient descent. At the M- th outer iteration, given the current transport matrix T (M) , we compute

where Cκ = f1(C g )pg1⊤ξ + 1η(pd )⊤f2(C d )⊤. When loss function is set as the square loss, f1(k) = k2 , f2(l) = l2 , h1(k) = k and h2(l) = 2l. Accordingly, the transport matrix at the (M+ 1) - th iteration is obtained by

It is easy to see that Eq.(9) can still be efficiently solved using the Sinkhorn-Knopp algorithm. Its closed-form solution is

where ¯K(M) = exp ( -L(Cg ,Cd , ¯ T (M) ) λ ), a and b are computed according to Eq.(5). The outer iterations are repeated until convergence or until the maximum number of iterations (set to 50) is reached, yielding the OT matrix ¯ T ∗ . Based on ¯ T ∗ , gene-disease pairs with the smallest transported mass are regarded as reliable negative pairs. To alleviate class imbalance, the number of selected negative pairs is set equal to that of positive pairs. The pseudo-code is given in Algorithm 2.

3.5 Node Representation Learning

According to the computation in Eq.(1), the nodes in graphs Gg and ˜Gd are respectively passed through two GCN encoders, and their node embeddings can be expressed as

Algorithm 2 Negative gene-disease pairs selection

Input:G(V g ,Eg ), ˜G(V d , ˜Ed ), pg, pd, λ, node features Xg ,Xd , max iterationsMmax.

Output: Negative gene-disease pairs E -.

• 1: Compute matrices C g ,C d and ¯ T (0) = pgpd T .
• 2: forM= 0 toMmax -1 do
• 3: Calculate L(C g ,C d , ¯ T (M) ) via Eq.(8).
• 4: L←L(C g ,C d , ¯ T (M) ).
• 5: ¯K ← exp(-L/λ).
• 6: Update optimal transport ¯ T (M+1) via Eq.(9).
• 7: end for
• 8: ¯ T ∗ ← ¯ T (Mmax) .
• 9: E - ← gene-disease pairs with smallest transport mass, |E -| = |E +|.

10: return E -

where Z g = [zg 1 , . . . , z g η ]⊤ and Z d = [zd 1 , . . . , z d ξ ]⊤, Xg and Xd denote the node feature matrices of genes and diseases, respectively. Ag and Ad are the corresponding adjacency matrices with self-loops. W g and W d are learnable parameters. We define the set of positive gene-disease pairs as ε+ = {(vg i , vd j )|Ahij = 1} and the set of negative pairs as ε- ⊂ {(vg i , vd j )|Ahij = 0}, Ah denotes the adjacency matrix of the gene-disease graph. The negative set ε- is constructed by selecting node pairs that receive the smallest transport mass in ¯ T, with the constraint |ε+| = |ε-|. Following the strategy in [Hadsell et al., 2006], the objective of node representation learning is to optimize W g and W d so that embeddings of connected node pairs are close in the embedding space, while embeddings of unconnected pairs are well separated. The objective function can be expressed as

where zg i and zd j are the embeddings of nodes vg i and vg j , respectively. The hyperparameter α balances the attraction of positive pairs and the repulsion of negative pairs, while γ > 0 is a margin that activates a penalty only when the distance between a negative pair falls below this threshold, thereby preventing representation collapse. The overall pipeline is shown in Figure 1.

3.6 Gene-Disease Associations Prediction

Here, the learned embeddings of gene nodes and disease nodes are denoted as ˜ Z g and ˜ Z d , respectively. For a given gene-disease node pair (vg i , vd j ), we obtain its joint representation by concatenating the corresponding node embeddings: ˜ zij = [˜ zg i , ˜ zd j ]. Each positive gene-disease pair is assigned a label y = 1, whereas each negative pair is assigned y = 0. We then employ a MLP as the classifier to distinguish between associated and non-associated gene-disease pairs (see Figure 2). Specifically, the concatenated representation ˜ zij is fed into the MLP to produce a prediction score: ˆ yij = σ(fMLP (˜ zij )). Let D denote the training set consisting of both positive and negative gene-disease pairs, with corresponding labels yij ∈ {0, 1}. The classification loss is defined as

4 Experiments and Results

4.1 Implementation Details

In GDAs-OT, we design two separate GCN encoders for the gene-gene graph and the disease-disease graph to learn node representations. Each encoder consists of a single GCN layer, with an input feature dimension of 768 and an output embedding dimension of 96. For GDAs prediction, we use an MLP composed of two Linear layers. The input dimension of the MLP is 96, and the output dimension is 2, corresponding to the binary classification task. The numbers of training epochs for the GCN encoders and the MLP are determined based on empirical analysis. The loss curves under different training epochs are shown in Supplementary Material Figure 3. According to the experimental results, the numbers of epochs are finally set to 20 for the GCN encoders and 15 for the MLP, respectively. All experiments were conducted on a server equipped with four NVIDIA A100 GPUs (80GB memory each) and an Intel Xeon Gold 6248R CPU.

4.2 Comparative Methods

To evaluate the effectiveness of GDAs-OT, we compare it with several recent approaches for GDAs prediction, including HOGCN, DGP-PGTN [Li et al., 2023], XGDAG, FusionGDA, GlaHGCL, and PhenoLinker [Andreu et al., 2025]. All comparative methods are implemented using the parameter settings reported in their original papers. Specifically, for comparison with XGDAG, we follow the original experimental configuration, use a subset of the DisGeNet dataset and select the same ten disease-related nodes as those used in the XGDAG experiments. This dataset contains a total of 3,208 gene and disease nodes and 4,302 association edges. In addition, since XGDAG formulates the task as a multi-class classification problem, we report Precision, Recall, and F1-score for the Positive class only to ensure a fair comparison. (see Table 3).

For the remaining comparison experiments, all methods were evaluated on the same dataset and under identical experimental settings. To maintain class balance, we randomly sampled 10,000 associated pairs and 10,000 non-associated pairs from the DisGeNet database as positive and negative samples for the binary classification task. The dataset was split into training and test sets with a ratio of 8:2, and all methods were assessed using the same evaluation metrics, including Accuracy, F1-score, Area Under the ROC Curve (AUC), and Area Under the Precision-Recall Curve (AUPR). All experiments randomly sampled samples from the DisGeNet dataset using 5 different seeds. The results are presented as the mean ± standard deviation. As shown in Table 4, GDAs-OT outperforms the competing methods across all evaluation metrics. In particular, Accuracy and F1-score are improved by more than 4%, while AUC and AUPR exhibit gains of nearly 2% and over 2%, respectively. These results indicate that GDAs-OT achieves superior overall predictive performance compared with existing approaches.

4.3 Ablation Experiment

To examine the contribution of different modules in GDAsOT, we conducted an ablation study. Specifically, we designed three model variants as follows. 1) The disease graph is used without structural expansion. 2) Based on the transport matrix ¯ T ∗ , select the gene-disease pair with the highest transport mass from the unassociated pairs as negative pairs. 3) Node representations for the gene-gene graph and diseasedisease graph are obtained directly from GCN encoders, without the subsequent node representation learning module.

These three variants are referred to as GDAs-OT(a), GDAs-OT(b), and GDAs-OT(c), respectively. The experimental results are reported in Table 5. Compared with the full model, all variants exhibit performance degradation to varying extents. Among them, GDAs-OT(b) shows a substantial drop across evaluation metrics, indicating that the proposed negative sample selection strategy based on the OT matrix plays a crucial role in model performance.

4.4 Parameter Sensitivity Analysis

In this section, we adopted 5-fold cross-validation on the training set to identify the optimal hyperparameter configuration. The final performance was reported as the average over the five runs. For all hyperparameters in the GDAs-OT model, we predefined candidate search ranges. The learning rate was selected from {0.0001, 0.001, 0.005, 0.007, 0.01}, α ∈ {0.1, 0.5, 1, 3, 5, 10}, β ∈ {1, 5, 7, 10, 50, 100} and the margin parameter γ ∈ {1, 5, 7, 10, 50, 100}. The crossvalidation results under different hyperparameter settings are shown in Supplementary Material Figure 4. The results indicate that excessively large or small values of lr, β and γ lead to a noticeable drop in accuracy. Interestingly, the model performance shows a negative correlation with the parameter α. A possible explanation is that overly large values of α place too much emphasis on separating negative pairs, which weakens the constraint that pulls positive pairs closer in the embedding space and ultimately degrades overall performance. Based on these observations, we set the final hyperparameters to lr = 0.001, α = 0.1, β = 5 and γ = 3.

4.5 Identify Genes Related to IrAEs

We selected three common types of irAEs from the existing research data [Ramos-Casals et al., 2020] for analysis: Psoriasis, Myasthenia gravis, and Guillain-Barre syndrome. All three irAEs were treated as diseases. We used all observed gene-disease pairs and an equal number of negative pairs to train GDAs-OT. After training, the model was applied to estimate association confidence scores between genes and each of the three irAEs. For each irAE, genes with confidence scores greater than 0.99 and ranked within the top-5 were selected as potential associated genes. The identified genes are listed in Table 6, where "Reported" indicates that the corresponding GDAs has been documented in the Dis- GeNet database. The results show that most genes identified by GDAs-OT have already been reported in DisGeNet, suggesting that the proposed method is effective in identifying disease-related genes. Notably, GDAs-OT also revealed novel GDAs absent from existing databases, including potential associations between psoriasis and NHSL1 and GFPT2, as well as between myasthenia gravis and SLC22A13 and PTBP3. Furthermore, recent studies have reported that GFPT2 is closely involved in immune-related functions such as immune cell infiltration and immune regulatory processes [Liu et al., 2024; Zhou and Dong, 2025].

5 Conclusion

In this work, we propose a GDAs prediction method based on optimal transport. By leveraging optimal transport, the proposed method identifies potential associations in the diseasedisease graph and reliable negative associations in the genedisease graph, thereby addressing a key limitation of existing graph-based methods that overlook false negative associations. Comparative experiments with several recent methods demonstrate that our approach consistently achieves superior performance across all evaluation metrics. In future work, we plan to further validate the irAEs-associated genes identified by our model through biological experiments, in order to assess their biological plausibility.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant No. 62462022, Grant No. 62572157).