CSGDN: Contrastive Signed Graph Diffusion Network for Predicting Crop Gene-phenotype Associations

We use datasets from three species including G. hirsutum you2023regulatory , B. napus tang2021genome and T. turgidum yang2024combined to input our model. For each species data, we perform TWAS process method to obtain associations between gene and phenotype. Initially, phenotype data is standardized using qqnorm function in R and Principal Component Analysis (PCA) is conducted for population structure data using TASSEL software bradbury2007tassel . FaST-LMM software lippert2011fast is employed to perform GWAS, considering phenotype data and population structure. We use Fusion software for performing TWAS gusev2016integrative to obtain TWAS $zscore$ of associations between partial genes and phenotypes. The positive or negative sign of TWAS $zscore$ indicates positive/negative regulation, respectively. Note that TWAS process can only calculate partial associations between genes and phenotypes. We separate the genes that can be associated with TWAS from those that cannot be calculated by TWAS within the species, and input them separately into the model.

G. hirsutum you2023regulatory . The phenotypes data include four phenotypes at the 4 DPA (days post anthesis) stage: Fiber length (FL), Fiber strength (FS), Fiber elongation rate (FE), Fiber Uniformity (FU). The reference genome is from Wang et al. wang2019reference . The result of TWAS process contains 523 associations between genes and FL phenotype, 1129 associations between genes and phenotype FS, 1521 associations between genes and phenotype FE, 509 associations between genes and phenotype FU.

B. napus tang2021genome . B. napus dataset includes 20 DPA and 40 DPA stages. phenotype type is only seed oil content (SOC). The reference genome (B. napus v4.1) can be downloaded from Genoscope (http://www.genoscope.cns.fr/brassicanapus/). The number of associations between genes in B. napus and SOC are 605 at the 20 DPA stage and 148 at the 40 DPA stage, separately. We regard the different stages as distinct phenotypes.

T. turgidum yang2024combined . T. turgidum dataset includes four phenotypes: Biomass (BM), Root length (RL), Root area (RA), and Root volume (RV). We use WEW v2.1 as the reference genome muslu2021comparative . The result of the TWAS process includes 36 associations with these four phenotypes, with BM having 1 association, RA having 2, RL having 17, and RV having 15.

As the complement to our model, we calculate the similarity between gene sequences as features. The software BLAST+/2.9.0 camacho2009blast+ is used for aligning sequences between gene pairs. BLAST can calculate the similarity scores between gene pairs. We use BLAST+ 2.9.0 camacho2009blast+ to obtain a similarity matrix $S\in\mathbb{R}^{n\times n}$, where $n$ = $|\mathcal{V}|$. We use the similarity matrix $S$ as the initial feature input for genes. To be specific, the input feature $h_{i}^{(0)}$ of gene $v_{i}$ is the $i$-th row of the similarity matrix.

The collection of large-scale crop sample data requires a high cost. Therefore, developing methods that can maintain high-precision prediction even with a low number of training samples is very necessary. In the field of biology, the effects of genes on phenotypes involves both positive and negative regulation, meaning that genes can either promote or inhibit phenotypes through complex mechanisms. This forms a bipartite signed graph, which consists of positive and negative edge types and two disjoint sets of nodes. In situations where crop sample data is scarce, the resulting signed graph structure is relatively sparse.

Traditional ranking models, such as PageRank page1999pagerank , can mine the potential relationships between nodes to increase the number of edges but are only suitable for graphs with a single type of positive edges. Wu et al. wu2016troll proposed Troll-Trust model (TR-TR) which is a variant of PageRank, without explanation of complex relationships between negative and positive edges. To better handle the signed graph with two edge types, we use a diffusion operation based on the Signed Random Walk with Restart (SRWR) algorithm proposed by previous researchers to uncover the potential associations between genes and phenotypes, obtaining a diffusion graph. The diffusion graph provides richer structural information, which can effectively alleviate the problem of sparsity in graph data, while also providing a new graph structure for subsequent contrastive learning operations.

To be specific, the diffusion method is based on the balance theory. In balance theory, there are four types of relationships between nodes: friend’s friend, friend’s enemy, enemy’s friend, and enemy ’s enemy. On this basis, we introduce a signed random surfer for bipartite signed graph. The surfer randomly surfs between nodes and traverses their relationships. When the surfer encounters a positive edge, it maintains a positive sign +; when it encounters a negative edge, it flips its sign to negative -. Initially, the surfer carries a positive sign + at node $u$. If the surfer is at node $v$ at this time and the restart probability is $c$, the probability of the surfer randomly surfing to neighboring nodes is $1-c$.

$\mathbf{r}_{u}^{+}$ is the probability that a positive sign surfer is at node $u$ from the seed node $s$, and $\mathbf{r}_{u}^{-}$ is the probability that a negative sign surfer is at node $u$ from the seed node $s$. We can formulate $\mathbf{r}_{u}^{+}$ and $\mathbf{r}_{u}^{-}$ by Equation 1 below where $\overleftarrow{\mathbf{N}}_{u}^{+}$ is the set of in-neighbors of node u , and $\overrightarrow{\mathbf{N}}_{v}^{+}$ is the set of out-neighbors of node v. Noted that node s is the seed node, and $\mathbf{1}(u=s)$ is 1 if u is the seed node s and 0 otherwise.

We use the adjacency matrix $\mathbf{A}$ of the signed graph $\mathcal{G}$ to vectorize Equation 1. Suppose $\mathbf{A}_{uv}$ is positive or negative for the signed edge $(u\rightarrow v)$ and 0 otherwise. $\mathbf{D}$ is the out-degree diagonal matrix, where $\mathbf{D}_{ii}=\sum_{j}|\mathbf{A}|_{ij}$. The semi-row normalized matrix is $\tilde{\mathbf{A}}=\mathbf{D}^{-1}\mathbf{A}$. We decompose $\tilde{\mathbf{A}}$ into two matrices: a positive semi-row normalized matrix ($\tilde{\mathbf{A}}_{+}$) and a negative semi-row normalized matrix ($\tilde{\mathbf{A}}_{-}$), such that $\tilde{\mathbf{A}}=\tilde{\mathbf{A}}_{+}-\tilde{\mathbf{A}}_{-}$.

Based on the above equations, we can vectorize Equation 1 as follows:

where $\tilde{\mathbf{A}}_{+}^{T}$ and $\tilde{\mathbf{A}}_{-}^{T}$ are the transposed adjacency matrices for positive and negative edges, respectively. $\mathbf{q}$ is the seed node vector at node $s$. Initially, set $\mathbf{r}^{+}=\mathbf{q}$, $\mathbf{r}^{-}=0$, and define $\mathbf{r}^{\prime}=\begin{bmatrix}\mathbf{r}^{+}\\ \mathbf{r}^{-}\end{bmatrix}^{T}$.

Since the simple balance theory with four types of relationships cannot explain the complex imbalanced relationships, introducing the balance attenuation factors $\beta$ and $\gamma$ is beneficial. When a negative walker encounters a negative edge at a node, its sign will switch to positive with probability $\beta$, or remain negative with probability $1-\beta$. Similarly, when the negative walker encounters a positive edge at node $m$, its sign will remain negative with probability $\gamma$, or switch to positive with probability $1-\gamma$. The diffusion model with balance attenuation factors becomes:

By repeated iterative computation of $\mathbf{r}^{+}$ and $\mathbf{r}^{-}$ , we can concatenate $\mathbf{r}^{+}$ and $\mathbf{r}^{-}$ into $\mathbf{r}=\begin{bmatrix}\mathbf{r}^{+}\\ \mathbf{r}^{-}\end{bmatrix}^{T}$. Then by computing the error $\delta$ between the current iteration result $\mathbf{r}$ and the previous iteration result $\mathbf{r}^{\prime}$, where $\delta=||\mathbf{r}-\mathbf{r}^{\prime}||$, we update $\mathbf{r}$ for the next iteration. The iteration stops when the error $\delta$ is smaller than a tolerance $\epsilon$.

For each node, the algorithm return an $\mathbf{r}^{+}$ and an $\mathbf{r}^{-}$. We combine all $r^{+}$ and $r^{-}$ into a positive matrix $\mathbf{r_{p}}\in\mathbb{R}^{n\times n}$ and a negative matrix $\mathbf{r_{n}}\in\mathbb{R}^{n\times n}$. Note that the relationship between edges ($u\to v$) and ($v\to u$) may differ, and the sign may even be opposite. This means $\mathbf{r_{p}}(u,v)\neq\mathbf{r_{p}}(v,u)$ may occur.

To address this and generate an undirected graph for subsequent process, we transpose $\mathbf{r_{p}}$ and $\mathbf{r_{n}}$ to $\mathbf{r_{p}}^{T}$ and $\mathbf{r_{n}}^{T}$ respectively where $\mathbf{r_{p}}^{T}(u,v)=\mathbf{r_{p}}(v,u)$ ( or $\mathbf{r_{n}}^{T}(u,v)=\mathbf{r_{n}}(v,u)$). For each node, we take the maximum value between $\mathbf{r_{p}}$ and $\mathbf{r_{p}}^{T}$ to form $\mathbf{r_{p}}_{\text{max}}$ and $\mathbf{r_{n}}_{\text{max}}$. Then, we calculate $\mathbf{r_{p}}_{\text{max}}-\mathbf{r_{n}}_{\text{max}}$ for each node at the corresponding position to generate the symmetric $\mathbf{r_{d}}$ matrix, which can also be called the diffuse matrix $\mathcal{S}$.

Generating different views is crucial in contrastive learning methods. In this study, we primarily focus on randomly removing edges on both the original graph and the diffusion graph to obtain different views. For each graph, we can construct a matrix $\mathcal{M}\in\mathbb{R}^{2\times n}$, where each column represents one edge (such as u->v) containing only existing edges. Then we generate a random masking matrix $\widetilde{\mathcal{R}}\in\mathbb{R}^{1\times n}$ drawn from a uniform distribution over $[0,1]$, denoted as $\widetilde{\mathcal{R}}_{i}\sim\text{Uniform}(0,1)$. Setting a threshold $p_{r}=0.1$, we reset $\widetilde{\mathcal{R}}_{i}=0$ to denote the deletion of the corresponding edge if $\widetilde{\mathcal{R}}_{i}<0.1$. The resulting matrix can be computed as: $\widetilde{\mathcal{M}}=\mathcal{M}\circ\widetilde{\mathcal{R}}$, where $(x\circ y)_{i}=x_{i}y_{i}$ is Hadamard product.

We perform this process twice on both the original graph $\mathcal{G}$ and the diffusion graph $\mathcal{S}$ separately, randomly removing 10% of the edges each time. Therefore, we can obtain four augmented graphs, denoted as $\mathcal{G}_{1}$, $\mathcal{G}_{2}$, $\mathcal{G}_{3}$, and $\mathcal{G}_{4}$, respectively.

After data augmentation, we obtain four augmented graphs including $\mathcal{G}_{1}$, $\mathcal{G}_{2}$, $\mathcal{G}_{3}$, and $\mathcal{G}_{4}$. For convenience, we define the set of augmented graphs as $\mathbf{G}_{k}=\left\{\mathcal{G}_{1},\mathcal{G}_{2},\mathcal{G}_{3},\mathcal{G}_{4}\right\}$, where k=1,2,3,4. Edge types are classified into positive and negative, meaning the two types of effects of genes on specific phenotypes, namely up-regulation and down-regulation. Consequently, it is imperative to design two distinct GNNs to separately aggregate information from positive and negative neighbors. We split each graph into two graphs containing only positive or negative edges, where $\mathbf{G}_{i}^{+}=(\mathcal{U}\cup\mathcal{V},\mathcal{E}^{+})$ and $\mathbf{G}_{i}^{-}=(\mathcal{U}\cup\mathcal{V},\mathcal{E}^{-})$. Utilizing a design akin to shu2021sgcl , node representations are learned from positive graphs using a positive GNN, and from negative graphs using a negative GNN. Parameters are shared in the same perspective for positive (negative) GNNs. GAT model velivckovic2017graph is used as the graph encoder and it computes as follows:

where $\zeta\in\left\{+,-\right\}$, $k$ represents the $k$-th augmented graph, $L$ denotes the number of GNN layers and $\mathcal{W}_{k}^{\zeta}$ is a learnable transformation matrix. $h_{i,k}^{(0),\zeta}$ denotes the input feature vector of the $i$-th node and $h_{i,k}^{(l),\zeta}$ is the representation of the $i$-th node for the $l$-th layer. $z_{i,k}^{\zeta}$ represents the final representation of $i$-th node in the $k$-th augmented graph. Note that we use two GAT layers here.

For those genes that lack TWAS associations, due to the lack of supervision from TWAS association information, the encoder trained specifically for genes in TWAS associations cannot be directly used for their encoding. To address this issue, we propose to train a multi-layer perceptron (MLP) to transform these genes that lack TWAS associations into a shared space with TWAS-associated genes:

where $h_{i}^{(0)}$ represents the input feature of gene $v_{i}$ in TWAS associations, and $z_{i}$ represents the final representation of gene $v_{i}$ after being encoded by the aforementioned TWAS framework. By minimizing the above Mean Squared Error (MSE) loss, we can learn the encoding capability of the main framework for TWAS-associated genes through MLP, which can be used for the encoding of genes lacking TWAS associations. Then we can obtain the representation of genes lacking TWAS associations through this MLP:

where $v_{i}\in\mathcal{V}_{\text{noTWAS}}$. The specific process is shown in Fig. 3.

Baselines. To validate the effectiveness of CSGDN, we compare our proposed model with several common methods in the fields of unsigned and signed graph neural networks.

$\bullet$ Unsigned GNNs. GCN kipf2016semi is a pioneering and notable GNN model tailored for unsigned graphs, featuring an effective layer-wise propagation mechanism. GAT velivckovic2017graph utilizes masked self-attentional layers, allowing nodes to attend to their neighbors’ features with varying weights without costly matrix operations or prior knowledge of the graph structure. GRACE zhu2020deep proposes a novel unsupervised graph representation learning framework that leverages contrastive objectives at the node level, creating two graph views through corruption, and maximizing the agreement of node representations in these views, utilizing diverse contexts via a hybrid scheme on structure and attribute levels, and demonstrating superior performance over state-of-the-art methods.

$\bullet$ Signed GNNs. SGCN derr2018signed utilizes balance theory to correctly aggregate and propagate the information across layers of a signed GCN model and generalizes GCN to signed graphs. SGCL shu2021sgcl introduces a novel graph contrastive representation learning techniques tailored for signed graphs, leveraging balance theory and dual contrastive strategies to achieve superior node representations across diverse datasets, including social and online gaming networks. SGNNMD zhang2022sgnnmd utilizes signed graph neural networks to predict deregulation types of miRNA-disease associations, achieving competitive performance by integrating structural and biological features from a signed bipartite network.

Hyper-parameters setting. we analyze the sensitivity of CSGDN to six key hyperparameters: $\alpha$, $\beta$, $feature$, $mask$, $predictor$, and $\tau$. The default configuration for each hyperparameter is as follows: $\alpha=0.8$, $\beta=0.01$, $feature=64$, $mask=0.4$, $predictor=\text{2-layer MLP}$, and $\tau=0.05$. These settings were determined based on the model’s overall highest AUC score. The AUC value is utilized as the primary metric to assess the sensitivity of the model to changes in these hyperparameters.

Reducing sample size. We randomly extract eighty percent of the G. hirsutum dataset. CSGDN shows the excellent results to face with small sample size.

Random Noise. To demonstrate that our model CSGDN has an outstanding performance when resisting interference, we achieve the effect of random perturbation to simulate noise by randomly flipping a certain proportion of edge signs in the G. hirsutum dataset. In this experiment, we set the proportion as 10% and 20%.

Task and evaluation metrics. We use AUC, Micro-F1,Binary-F1 and Macro-F1 to evaluate the results on the link sign prediction task. For each dataset, we randomly split edges into a training set and a testing set with a ratio 8:2. Note that superior performance is indicated by higher values for all these four evaluation metrics.

Performance of CSGDN compared with baselines (Q1). To answer Q1, we compare CSGDN with current state-of-the-art methods. We primarily use two types of GNN frameworks as baselines: unsigned GNN and signed GNN. For unsigned GNNs, we used GCN, GAT, and GRACE, while for signed GNNs, we employed SGCN, SGCL, and SGNNMD. We use link prediction as the evaluation task, with AUC, Binary-F1, Micro-F1, and Macro-F1 as the evaluation metrics for model performance. Across three common crop datasets (G. hirsutum, B. napus, and T. turgidum), the evaluation metrics of CSGDN generally outperform those of the state-of-the-art baselines. As shown in Table 2, CSGDN demonstrates strong performance in the link prediction task on crop datasets.

Performance of CSGDN when addressing small sample size and random noise (Q2). As shown in Tables 3 and 4, we address Q2 by verifying that CSGDN can effectively overcome two issues: the costs of samples and noise. For the first issues we randomly reduce the sample size of the G. hirsutum dataset to 80% and subsequently divide this dataset into training and testing sets. The results presented in the Table 3 indicate that CSGDN outperforms all baselines on G. hirsutum datasets with randomly reduced sample sizes, demonstrating its effectiveness in handling small sample datasets. This suggests that we can reduce experimental costs and durations by minimizing the sample size while still achieving excellent prediction outcomes. Then for noise, we utilize two G. hirsutum datasets with the proportion of 10% and 20% perturbations. As shown in Table 4, our model CSGDN outperforms the vast majority of baselines. This reflects that our model has strong anti-interference ability against various types of noise through the contrastive learning methods.

We conduct ablation study to assess the effectiveness of different components in our proposed model to answer Q3. In this subsection, we employ sign perturbation as the graph augmentation method to analyze performance. Specifically, we compare CSGDN with its three variants: $\text{CSGDN}_{\mathrm{w/o\ diffuse}}$, $\text{CSGDN}_{\mathrm{w/o\ aug}}$, and $\text{CSGDN}_{\mathrm{w/o\mathcal{L}_{CL}}}$, which are defined as follows:

• •

  $\text{CSGDN}_{\mathrm{w/o\ diffuse}}$ : The graph diffusion in contrastive learning is removed, and instead, we utilize the same two original graphs without diffusion step for the next augmentation step.

• •

  $\text{CSGDN}_{\mathrm{w/o\ aug}}$ : The graph augmentation step is removed. In this variant, graphs $\mathcal{G}$ and $\mathcal{S}$ instead of augmented graphs $\mathbf{G_{K}}$ are exploited during training.

• •

  $\text{CSGDN}_{\mathrm{w/o\mathcal{L}_{CL}}}$ : The part of contrastive learning is removed and ignores the contrastive loss.

As shown in Table 5, we demonstrate that all three components are essential to CSGDN’s performance. Each components plays a unique role in enhancing the model’s ability.

To answer Q4, we analyze the sensitivity of CSGDN to six key hyperparameters: $\alpha$, $\beta$, $feature$, $mask$, $predictor$, and $\tau$.

The hyperparameter $\alpha$, shown in Fig. 4(a), serves as the weight coefficient that balances the significance of the inter-view contrastive loss and intra-view contrastive loss. We evaluate the model’s performance with $\alpha$ set to values from the set ${0.2,0.4,0.6,0.8,1.0}$. We find that the model performs better when $\alpha$ is between 0.5 and 0.8, while anything outside of this range decreases performance. Therefore, for this dataset intra-view, i.e., learning a consistent representation across augmented graphs, is more important.

The hyperparameter $\beta$, shown in Fig. 4(b), controls the trade-off in the joint loss function between the link sign prediction loss ($\mathcal{L}_{label}$) and the contrastive learning loss ($\mathcal{L}_{CL}$). We vary $\beta$ over the set ${0,0.1,0.01,0.001,0.0001,0.00001}$ and observe that the highest AUC score is achieved when $\beta=0.01$. We find that the model performance tends to increase when $\beta>0.01$, and the model reaches its best performance at $\beta=0.01$, followed by a decrease in performance. This illustrates the importance of contrast learning in the model, which decreases when the contrast learning effect on model performance is extreme, i.e., when the $\beta$ value is too small or too large.

The hyperparameter $feature$, shown in Fig. 4(c), refers to the node embedding dimension. We test six different node embedding dimensions, ranging from 8 to 128, and evaluate their impact on model performance. We find that the model performance decreases extremely when $feature=128$, which may be due to feature sparsity caused by the elevated dimension of the feature space, which prevents the model from effectively establishing correlations between nodes.

The hyperparameter $mask$, shown in Fig. 4(d), represents the ratio of edges randomly dropped from both the original and diffusion graphs. We explore the values of $mask$ in the set ${0,0.2,0.4,0.6,0.8}$ and find that the model performs optimally when the mask ratio is set to 0.4. We find that the model performance fluctuates as mask ratio increases, but generally performs better when mask ratio is small. When mask ratio is large, too much information is lost, resulting in incomplete information and lower model performance.

The hyperparameter $predictor$, shown in Fig. 4(e), determines the best architechture of the downstream link prediction task, allowing for variations among a 1-layer MLP, 2-layer MLP, 3-layer MLP, and 4-layer MLP. The number of layers directly influences the model’s prediction accuracy. We found that as the number of layers of $predictor$ increases, the model performance first increases and then decreases, and reaches the highest value at 3-layer MLP, while 4-layer MLP may be due to the overfitting of the model because of the excessive number of layers of the neural network, which is ultimately manifested in the decrease of the model performance on the test set.

Finally, the hyperparameter $\tau$, shown in Fig. 4(f), is the temperature parameter used in the contrastive loss function to regulate the similarity contrast between positive and negative edges. We test $\tau$ values from the set ${0.05,0.1,0.2,0.4,0.8}$, and the default value of $\tau=0.05$ yields the best performance. We find that the model performance is better when the $\tau$ value is lower, and generally lower when $\tau>0.05$, which also indicates that this dataset needs to learn a consistent representation as much as possible in the comparative learning, and smoothing operation here will reduce the performance.

In this section, we initiate a case study focusing on G. hirsutum genes that associated with the four types of phenotypes including Fiber Elongation rate (FE), Fiber Uniformity (FU), Fiber Strength (FS), Fiber Length (FL). By leveraging the optimal hyperparameter configurations outlined, we utilize the associations between G. hirsutum genes and four types phenotypes predicted by TWAS process as the training set to train the CSGDN’s TWAS Frame. Then, the randomly-selected gene-phenotype associations types that TWAS cannot calculate in G. hirsutum are used as irrelevant associations. CSGDN can predict probabilities of three types associations, including up-regulation, down-regulation and irrelevant associations. We take the type of the highest predicted probability as the current gene and phenotype type. Thus, for the target four phenotypes, we can predict association types with target G. hirsutum genes. As shown in Table 6, we randomly list 10 associations of each phenotype. We use the result associations from TWAS as the reference for our predictions. For example, for phenotype FE, 2 out of 16 predictions are incorrectly compared with TWAS results. For phenotype FU, 2 out of 9 predictions are incorrect from the TWAS results. Therefore, the case studies demonstrate the usefulness of CSGDN in discovering novel gene-phenotype associations and can be validated by TWAS results.