drGAT: Attention-Guided Gene Assessment of Drug Response Utilizing a Drug-Cell-Gene Heterogeneous Network

1. Introduction

Drug discovery is an expensive and lengthy process with many obstacles [1,2]. The difficulty of this process arises from the efforts needed to ensure a compound is both safe and effective. Sensitivity to a drug involves the mechanism of the drug compound and a complex interplay of various factors internal and external to a cell. These factors include the cellular context (or state) that is determined by the repertoire of transcripts and proteins, alterations to this repertoire due to disease, and the interactions between these components [3,4]. Biomarkers are crucial in drug development to understand the best use of a compound and aid our understanding of the disease biology [5,6,7]. Machine learning (ML) has emerged as an approach to understanding the role of particular genes in drug response. This makes a huge impact on the efficiency and success rates of developing new drugs [8]. The utilization of ML in this domain leverages large amounts of biological and chemical data, helping to identify drug candidates and predict their efficacy efficiently. When it comes to understanding biological phenomena, ML algorithms, especially deep learning (DL) models, have demonstrated their capability to predict the complex patterns and relationships within biological data [9,10].

However, despite the advances in using machine learning for drug discovery, a critical difficulty remains; which is the “black box” nature of the methods. While DL models have made huge strides in identifying patterns and predicting outcomes within several areas, their internal decision-making process is unclear, leading to concerns about their trustworthiness and reliability.

Interpretability in machine learning has received significant recent consideration [11,12,13] including in the area of drug discovery [14,15]. The attention mechanism, introduced in the Transformer architecture [16], is being widely used due to its capacity to include “attention coefficients” as trainable parameters to optimize how much attention the model assigns to individual components, such as a word in a sentence (i.e., a linear sequence of data); these attention parameters can be used to interpret the model’s output. While words in a sentence express a single coherent idea, regulatory networks biology are more dependent on interconnections between components; this makes graph data structures a useful representation in biology. Graph Neural Networks (GNNs) have been developed for handling this type of data. GNNs aggregate node-level features in local graph neighborhoods, thereby forming features that also incorporate node relations. GNNs that make use of the attention mechanism are known as Graph Attention Network (GAT) [17]. Thus, GAT leverages the attention mechanism to generate attention weights, which weigh the influence of feature representations from neighboring nodes when deriving node-level features. By utilizing GAT, we can offer insights into the importance of genes and enhance the interpretability of drug response models by analyzing neighborhood relationships.

We propose a novel interpretable deep learning method called “drGAT”, which leverages GAT to process a large-scale heterogeneous network. Our model processes a heterogeneous graph composed of drug compounds, cell lines, and genes. This heterogeneous graph has been constructed using integrated drug screening and molecular profiling data from the NCI60 cell line pharmacogenomics dataset (drug screening data collected on 60 cancer cell lines) taken from CellMinerCDB [18] via rcellminer [19]. The GAT layer is employed to learn embedding representations from the interconnected structure within this heterogeneous network. This approach gives us knowledge of the relative significance of individual components [20,21] and their contributions to the overall model.

We utilize drGAT for two tasks: i) drug-cell line association prediction and ii) interpretability of the individual gene importance to prediction. We make use of multi-task learning, whereby multiple learning tasks are solved simultaneously to benefit from the shared patterns and distinctions among these tasks. Towards this aim, this study focuses on drug compounds with a DNA-damaging mechanism of action (269 drugs) available in the NCI60 of which there are more than in other comparable datasets. Regarding drug response prediction, drGAT has demonstrated superior performance over existing models. Additionally, this model also shows high prediction accuracy on data not present in NCI60 but found within the GDSC dataset.

In addition, the attention coefficients from drGAT are utilized for interpreting the predictor (e.g., gene) importance and relation among components of the heterogeneous network nodes (i.e., drugs, cell lines, and genes). We present an analysis of attention coefficients for the subset of nodes that correspond to genes that propose explanations for the response mechanisms of particular drugs. Moreover, relevant biological processes for analyzed compounds are also described using the attention coefficients and an over-representation analysis.

In summary, our model shows high accuracy for drug response prediction, allowing us to determine the predictor importance using attention coefficients. This interpretability allows users to explore the relationship among drug structures, cell lines, and genes, providing insights for further drug development investigation.

The main contributions include:

• We create a heterogeneous graph involving DNA-targeting drugs that includes three entity types: genes (i.e., gene expression), cell lines (i.e., drug responses), and drugs (i.e., structures) using the NCI60 [22].
• We evaluate the drug-gene associations suggested by attention coefficients and systematically compare these to existing scientific publication text. We examine the abstracts of journal papers for co-mentions of drug-target relationships from our results.
• Our results propose drug-cell line sensitivity associations based on the attention coefficients derived from the GAT. These associations were then validated through comparison with independent experimental data.

2. Results

2.1. Model Performance

In Table 1, we present the results of the drug response prediction. A comparison is made between four baseline methods and our model with several GNN layers. We train and fit the model 5 times and values show the average and standard deviation. Notably, GATv2 achieved an accuracy of 0.779, a precision of 0.775, a recall of 0.741, and an F1 score of 0.757, outperforming other methods in terms of accuracy and F1 score, two of the most important metrics in binary classification. We utilize GATv2 as the GNN layer in our model for further investigation. Furthermore, our model can provide interpretations based on attention coefficients that existing models cannot achieve.

2.1.1. Analysis of GNN Architectures and Their Impact on Performance

Three key differences exist among GNN architectures: convolution, attention, and long-range dependencies. These differences significantly impact each model’s performance and capabilities.

The MPNN (Message Passing Neural Network) updates node features by receiving messages from neighborhood nodes, capturing complex relationships within the graph. However, its computational intensity and risk of over-smoothing with increased layers limit its performance, with an accuracy of 0.620 ± 0.126.

GCN (Graph Convolutional Network) uses graph convolution operations to update node features, ensuring low computational cost. However, GCN averages neighboring nodes’ information, which can result in losing important features. This limits GCN’s performance to an accuracy of 0.678 ± 0.036.

GAT (Graph Attention Network) uses attention mechanisms to weigh the importance of neighboring nodes, focusing on relevant features. This enhances GAT’s interpretability and performance, with an accuracy of 0.763 ± 0.018, though it increases computational cost.

GATv2 is an enhanced version of GAT, with improved attention mechanisms for more accurate feature learning and better scalability. Despite sharing GAT’s higher computational cost, it achieves the highest performance among the models, with an accuracy of 0.779 ± 0.002.

The Graph Transformer applies the transformer model to graph data, capturing long-range dependencies and handling complex graph structures. However, its high computational cost limits its performance to an accuracy of 0.764 ± 0.01, not surpassing the efficiency and effectiveness of GATv2.

2.2. Interpretation

To interpret the results further, we extracted the attention coefficients from two GAT layers. Each layer includes 5 attention heads. For ease of use, we simplified the shape of the attention matrices from (n, n, 5) to (n, n) by averaging them. We then added these matrices together and used the resulting matrix in subsequent steps.

2.2.1. Evaluation of Drug-Gene Associations in Scientific Literature

We conducted an analysis of our drug-gene co-occurrences in PubMed abstracts using, for each drug, the top 5 genes by attention coefficient magnitude. We utilized the NCBI ESearch tool to retrieve abstracts that mention both the drug and gene of each relationship [23]. Queries used the following format: https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=pubmed&term="DRUG"&term="GENE".

Figure 2 displays a heatmap illustrating the drug-gene co-occurrences based on PubMed abstracts. Out of 1,345 drug-gene attention-based relationships, 370 of the drug-gene co-occurrences were found in the abstracts from 10,100 articles. The color corresponds to the number of related publications by log scale. The symbols in Figure 2 indicate different relationship types. The ♡ symbol represents relationships predicted by drGAT and present in the DTI dataset (53 instances). The ♣ symbol indicates relationships predicted by drGAT but not present in the DTI dataset (495 instances). The ♢ symbol denotes relationships only present in the DTI dataset (2 instances). Our model identified 53 co-occurrences in 55 (96.4 %), demonstrating high coverage of known DTIs in the dataset. Regarding predicted-only DTIs (♣), 156 out of 495 total predictions (31.5 %) have at least one PubMed abstract. This indicates our model’s potential to suggest new drug-gene relationships.

Figure 1. drGAT Overview. A heterogeneous graph is constructed using drug structures, drug response, gene expression data. Each input type is transformed using RBF kernel followed by three linear layers to standardize their feature sizes, making them compatible for integration. These transformed inputs are then combined into a unified feature matrix, which is used as input to a GAT layer with an adjacency matrix. For prediction, the GAT layer output is concatenated to align with the drug response data. Following a fully connected layer, a sigmoid function generates predicted values $\widehat{y}$, which are then inputs for the binary cross-entropy loss. To assess interpretability, attention coefficients are used with over-representation analysis and predicting untested sensitivities.

Figure 1. drGAT Overview. A heterogeneous graph is constructed using drug structures, drug response, gene expression data. Each input type is transformed using RBF kernel followed by three linear layers to standardize their feature sizes, making them compatible for integration. These transformed inputs are then combined into a unified feature matrix, which is used as input to a GAT layer with an adjacency matrix. For prediction, the GAT layer output is concatenated to align with the drug response data. Following a fully connected layer, a sigmoid function generates predicted values $\widehat{y}$, which are then inputs for the binary cross-entropy loss. To assess interpretability, attention coefficients are used with over-representation analysis and predicting untested sensitivities.

2.2.2. Evaluation of Predicted Drug-Cell Line Response by Comparison with GDSC

We systematically evaluated our results by comparing them with the Genomics of Drug Sensitivity in Cancer (GDSC) dataset [24]. The GDSC dataset is an extensive collection of drug response screening, which includes 978 cell lines and 542 drug compounds. The NCI60 subset we use here and GDSC overlap with 55 cell lines and 201 drugs. For evaluation, we use 44 overlapping drugs that have GDSC but no NCI60 response data.

To assess these drug-cell line responses, we focused on their attention coefficients using our attention model, ensuring an accurate interpretation of their relationships. Our evaluation involved defining drug-cell sensitivity associations based on the attention coefficients of genes. We computed the cosine similarity between the attention directed towards genes from drugs and cell lines, selecting the cosine similarity due to its effectiveness in comparing sparse vectors. A cosine similarity value exceeding 0.5 was classified as positive (indicating sensitivity), whereas values below this threshold were deemed negative (resistant) (Figure 3A).

The GDSC dataset provides IC50 in micromolar ($\mu M$) units. IC50 is described as the "half-maximal inhibitory concentration," referring to the concentration of a substance (such as a drug or inhibitor) required to inhibit viability or a specific biological function by 50%. This enables us to determine drug sensitivity based on dosage. We calculated four key metrics: accuracy, precision, recall, and the F1 score, across various dose thresholds ranging from 1 to 1400 $\mu M$ (Figure 3B). Figure 3B illustrates the impact of varying the sensitivity threshold value on the performance metrics of a classification model: accuracy, precision, recall, and F1 score. When the threshold is set low, the model achieves high accuracy and precision, identifying positive cases and minimizing false positives. However, a low threshold also leads to higher false negatives, as indicated by the lower recall. This scenario suggests that the model is conservative in predicting positive outcomes, leading to many positives being classified as negatives.

Conversely, setting a higher threshold improves recall, meaning the model is better at identifying all positive cases, including those that are harder to classify correctly. This comes at the expense of accuracy and precision as the model becomes more lenient, increasing the likelihood of false positives.

The balance between these metrics—accuracy, precision, and recall—highlights the trade-off in setting the threshold. While each drug may be effective at different concentrations, lower IC50 concentrations tend to be more clinically relevant. The threshold where each metric (accuracy, precision, recall, and F1 score) achieves a balance for this analysis is around 30 $\mu M$. (Figure 3B). At this threshold, the metrics are as follows: Accuracy: 0.62, Precision: 0.70, Recall: 0.8, and F1 Score: 0.75. This optimal threshold is one that balances performance across all four metrics, managing the trade-offs between identifying positive cases, minimizing false positives, and accurately classifying negative cases. Lower IC50 values can reflect the clinical relevance for compounds utilized in the analysis. The average IC50 value from the GDSC dataset is 100$\mu M$. The discrepancy between the average IC50 value and the IC50 threshold used for analysis may arise because of the binarization to distinguish between positive and negative cases.

Such insights are crucial for tuning classification models to achieve the desired balance between sensitivity (recall) and specificity (accuracy and precision), especially in applications where the cost of false positives and false negatives varies significantly.

2.2.3. Drug-Target Interactions Assessment from Attention Coefficients

The intricate network of drug-target interactions and quantifying their relationships from the model results through attention coefficients can aid our understanding of mechanisms of action. Figure 4 (A) presents a visualization of these drug-gene relationships, employing the drGAT model to highlight the interplay between chemotherapeutic drugs and their genetic targets.

Doxorubicin (NSC-759155 and NSC-123127) and daunorubicin (NSC-756717) are both anthracycline chemotherapeutic agents and are known as TOP2 inhibitors [25,26]. Results from drGAT retain this known drug-target association; TOP2A has the highest attention coefficient for both drugs. Doxorubicin and daunorubicin share an association with KRT14; KRT14 is a gene that is a member of the keratin 1 family. There is evidence that keratin expression is altered in response to doxorubicin treatment [27] and are altered in doxorubicin-resistant cells (including specifically, KRT14) [28]. Additionally, there is interest in targeting keratins through peptide-drug conjugate featuring doxorubicin [29]. Separately, SERPINA6, returned as a top 5 attention coefficient gene for doxorubicin, has been utilized as a resistance marker gene for doxorubicin [30].

Topotecan (NSC-759263), Camptothecin (NSC-94600), and Irinotecan (NSC-759878) are topoisomerase inhibitors, specifically targeting the TOP1 cleavage complex [7,31,32], and the TOP1 drug-gene association appears in the drGAT results. MMP3 was returned as a drug-gene association for topotecan by drGAT. MMP3, and other matrix metalloproteinases (MMPs), have been analyzed as a part of a Phase I dosing study of intraventricular topotecan on children with neoplastic meningitis due to reports of correlation with leptomeningeal disease [33].

While not a comprehensive analysis, the above results show that our model retains high attention coefficients of known drug-target interactions fed as training data, and also, predicts several known drug-gene relationships that are of research interest. This suggests that our predicted drug-gene associations can be indicative of potential targets or drug-gene relationships appropriate for further study.

2.2.4. Over-Representation Analysis with Attention Coefficients

We worked to elucidate the functional roles of genes associated with the various drugs, leveraging the results provided by the attention coefficients derived from our model. We conducted an over-representation analysis (ORA) using the attention coefficients attributed to the genes associated with the drugs. Specifically, genes possessing attention values greater than 0 for each drug were taken as input for ORA. We conducted the ORA analysis using gseapy [34] and the MSigDB Hallmark 2020 [35] collection of 50 gene sets accessed via gseapy. We highlight the biological processes with p-values adjusted by Benjamini-Hochberg less than 0.05 in Figure 4B, which depicts associations with different hallmark biological processes, with the cumulative count of drugs linked to each function.

Our study focused on drugs related to DNA damage. Consequently, our analysis reproduced processes represented in the Hallmark gene sets commonly associated with these drugs such as response to ultraviolet (UV) radiation, E2F targets (i.e., cell cycle-related targets of E2F transcription factors), xenobiotic metabolism, and DNA repair. The E2F targets are crucial in cell cycle regulation and may represent how these drugs impact proliferation [36]. Similarly, the involvement of xenobiotic metabolism pathways highlights how these drugs are processed. Related, ABCB1 was a recurrent overalapping gene in the "KRAS Signaling Up" gene set; this drug efflux protein is responsible for reducing concentrations of toxic compounds, such as cancer medications, and it has been shown to be associated with the resistance to topoisomerase inhibition [37,38,39]. The gene set "Allograft Rejection" was identified with the largest number of compounds. While not obvious from the name "Allograft Rejection" (i.e., transplant reject), this gene set contains many key cancer-related genes, including various kinases (e.g., EGFR, ITK, MAP4K1, LCK). Of these genes, one recurrent overlapping gene was EGFR for 80 of topoisomerase inhibitors (both TOP1 and TOP2; 130 topoisomerase inhibitors total); the various modes of interactions between topoisomerases and EGFR have been previously reported [40]. These results corroborate our expectations of the molecular mechanisms affecting response to these drugs.

3. Methods

3.1. Building Input Matrix

In this section, we explain the construction of the heterogeneous graph. The resulting graph incorporates drug structures, cell line drug response, gene expression data, and their known drug-target interactions and similarities.

3.1.1. Preprocessing

First, we selected a subset of data from NCI60 [22] via rcellminer [19]. This dataset comprises gene expression, drug response, and drug structure data.

The complete drug response matrix contained 23,191 drugs for 60 cells. From this, we selected 269 drugs with a mechanism of action related to DNA damage [41]. Specifically, the drug response matrix forms $X_{DR}\in {\mathbb{R}}^{n\times m}$, where n is the number of drugs and m is the number of cells (i.e., $n=269$, $m=60$).

The full NCI60 gene expression data includes 23,826 gene transcripts, where we selected 2,383 genes (the top 10% of genes with the greatest standard deviations). The selection of this highly variable gene set should allow better discrimination between different observations in the dataset. We also merged another group of genes (281 genes) that are involved in drug-target interactions. This combination resulted in a unique set of 2,718 genes. The gene expression data is part of matrix $X_{GE}\in {\mathbb{R}}^{m\times l}$, where $m(=60)$ and $l(=2,718)$ denote the number of cell lines and genes, respectively.

Subsequently, we collected the drug-target interaction data from DrugBank [42]. This dataset has 19,017 drug-target interactions (DTI) with 7,756 unique drugs and 4,755 genes. We selected 100 unique National Service Center (NSC) numbers (NSCs represent batches of screening with a particular drug structure) and 403 genes, with drugs and genes overlapping with the NCI60 pharmacogenomic data. To be the same size as other data, the matrix $X_{DTI}$ is filled by 0 for the unknown DTI. Therefore, $X_{DTI}$$\in {\mathbb{R}}^{n\times l}$ (i.e., $n=269$, $l=2,718$).

In addition, we also created the drug fingerprints representation matrix. SMILES structures are converted to the Morgan fingerprints [43] using RDKit [44]. We create a vector of 2048 length for each drug and concatenate to make a matrix, denoted as $X_{MF}\in {\mathbb{R}}^{n\times o}$, wherein o represents the number of drug features (i.e., $n=269$, $o=2,048$).

3.1.2. Feature Matrix

To make a feature matrix, we utilize similarity matrices as inputs created from $X_{DR},X_{GE},$ and $X_{MF}$. Preceding the construction of the heterogeneous graph, similarity matrices were individually constructed for cell lines ($S_c$), genes ($S_g$), and drugs ($S_d$). These matrices capture the element-wise similarity to homogenous entities such as gene-gene, cell-cell, and drug-drug similarity relationships.

The RBF Kernel [45] was employed to create a similarity matrix, defined as follows.

$$\begin{array}{cc}\hfill S_{ij}& =exp\left(-\gamma \left|\right|X_i-X_j{\left|\right|}^2\right),\hfill \end{array}$$

(1)

where S is a similarity matrix, $\gamma$ denotes a hyperparameter (we utilize ${\displaystyle \frac{1}{\mathrm{length}\left(X_i\right)}}$), i and j are indexes, and X is an input matrix.

From Equation (1) with $X_{DR}$, we obtained the cell similarity matrix $S_c$, and the gene similarity matrix $S_g$ was built using Equation (1) based on the gene expression data $X_{GE}$. Then, Equation (1) was applied to obtain the similarity matrix to the $X_{MF}$ to create a drug similarity matrix $S_d$.

As a result, three similarity matrices were obtained: $S_c\in {\mathbb{R}}^{m\times m}$ for cell lines, $S_d\in {\mathbb{R}}^{n\times n}$ for drugs, and $S_g\in {\mathbb{R}}^{l\times l}$ for genes. These matrices have different dimensions. To address this, we utilize three linear layers to transform them to the same size of hidden units. The unified feature matrix X is defined as follows:

$$X=\left[\begin{array}{ccc}{\mathrm{Linear}}_d\left(S_d\right)& 0& 0\\ 0& {\mathrm{Linear}}_c\left(S_c\right)& 0\\ 0& 0& {\mathrm{Linear}}_g\left(S_g\right)\end{array}\right],$$

(2)

where each linear layer has a different input size and the same output size. The output size is a hyperparameter.

3.1.3. Adjacency Matrix

For the input matrix, we also create the adjacency matrix from $X_{DR},X_{GE},$ and $X_{DTI}$. The drug-cell association matrix $A_{dc}$ is derived from the drug response matrix, $X_{DR}\in {\mathbb{R}}^{n\times m}$. Values in this matrix have already been normalized to z-scores with an average value of 0 and a standard deviation of 1. Then, we applied the following procedures to make it a graph structure. If the z-score in the drug response data for the NCI60 is greater than 0, we consider that the cell line is sensitive to the drug, and the value is 1. Otherwise, the cell line is resistant, and the value is set to 0.

Subsequently, we establish a cell-gene association matrix, $A_{cg}$, derived from the gene expression dataset $X_{GE}\in {\mathbb{R}}^{m\times l}$. Our initial step includes standardizing the data across columns (cell lines). After this, the cell-gene association matrix $A_{cg}$ is generated from $X_{GE}$, with processing similar to drug response matrix $X_{DR}$.

Lastly, we create the association matrix between drugs and genes from the drug-target interaction matrices $X_{DTI}$. To make a matrix with binary values, datasets with numerical values are binarized with values greater than 0 set to 1; otherwise, they are set to 0. Then, all datasets are merged as the drug-gene association matrix $A_{dg}\in {\mathbb{R}}^{n\times l}$ whereby if an interaction is present in any of the datasets, it will be retained in the final merged association matrix.

Based on the three matrices, a unified adjacency matrix A is created as the following Equation (3).

$$A=\left[\begin{array}{ccc}0& A_{dc}& A_{dg}\\ A_{dc}^{⊺}& 0& A_{cg}\\ A_{dg}^{⊺}& A_{cg}^{⊺}& 0\end{array}\right].$$

(3)

$A\in {\mathbb{R}}^{(n+m+l)\times (n+m+l)}$ represents the combined drug, cell, and gene components, where the dimension of F is (4430, 4430). In this process, the adjacency matrix A with non-zero values is represented as 1; otherwise, 0. From these processes, the input matrix X and the adjacency matrix A are utilized for the input of the GNN model.

3.1.4. Preventing Data Leakage via Masking

For evaluation, the drug response matrix $X_{DR}$ is partitioned into distinct subsets: training, validation, and test data (for precise specifications, refer to Appendix B.1). The test dataset is a subset of the drug-cell association matrix $A_{dc}$. This overlap raises concerns about potential data leakage. To counteract this issue, we mask the association values linked with the test data. Specifically, we set the values of test data combinations to 0 in the $A_{dc}$ matrix. This masking process affects 20 % of total combinations, which corresponds to 3,228 drug-cell line pairs. We implement this masking utilizing the following equation:

$${\left(A_{dc}\right)}_{ij}=\left\{\begin{array}{cc}0,\hfill & \mathrm{if}(i,j)\mathrm{is}\mathrm{a}\mathrm{test}\mathrm{data}\mathrm{index},\hfill \\ {\left(A_{dc}\right)}_{ij},\hfill & \mathrm{otherwise},\hfill \end{array}\right.$$

(4)

where i is the index of a drug, and j is the index of the cell for test data.

3.2. drGAT Model

This model consists of two $\mathrm{GAT}$ layers and a single fully connected layer. For the implementation of our drGAT model, we utilized Graph Attention Network v2 (GATv2) [46] as the core graph neural network layer. GATv2 was chosen due to its improved attention mechanisms, which allow for more expressive feature transformations compared to the original GAT. This choice was made to enhance the model’s ability to capture complex relationships within our heterogeneous graph structure.

Equation (5) illustrates the first block of the GAT layer with feature matrix X and adjacency matrix A. The subsequent steps involve adapting the model’s output to predict drug sensitivity in a binary manner.

$$Z_1=\mathrm{Dropout}\left(\mathrm{ReLU}\left(\mathrm{GraphNorm}\left(\mathrm{GAT}(X,A)\right)\right)\right),$$

(5)

where $\mathrm{Dropout}$ denotes the dropout layer, $\mathrm{ReLU}$ is a ReLU activation function, $\mathrm{GraphNorm}$ is a graph normalization layer [47], and $\mathrm{GAT}$ is a Graph Attention Layer.

Here, we utilized Graph Attention Layer (GAT) [17] defined as follows:

$$\mathrm{GAT}\left(x_i\right)={\alpha }_{ii}{W}_s{x}_i+\sum _{j\in \mathcal{N}\left(i\right)}{\alpha }_{ij}{W}_t{x}_j,$$

(6)

where $x_i$ is an input, $\alpha$ is the attention coefficient matrix, $W_s$ and $W_t$ are weight matrices, and $\mathcal{N}\left(i\right)$ is some neighborhood of node i in the graph. Here, the attention coefficients ${\alpha }_{ij}$ are computed as

$${\alpha }_{ij}={\displaystyle \frac{exp\left(LReLU\left(a_s^{⊺}{W}_s{x}_i+a_t^{⊺}{W}_t{x}_j\right)\right)}{{\sum }_{k\in \mathcal{N}\left(i\right)}exp\left(LReLU\left(a_s^{⊺}{W}_s{x}_i+a_t^{⊺}{W}_t{x}_k\right)\right)}},$$

(7)

where $a_s$ and $a_t$ are weight vectors, and LReLU describes LeakyReLU. This attention coefficient matrix $\alpha$ is utilized for further interpretability assessment.

From $Z_1$, we utilize the same structure and obtain $Z_2$ as follows.

$$Z_2=\mathrm{Dropout}\left(\mathrm{ReLU}\left(\mathrm{GraphNorm}\left(\mathrm{GAT}(Z_1,A)\right)\right)\right).$$

(8)

Drug-cell line associations are predicted by concatenating their respective embeddings (from $Z_2$ ), and passing this input through a fully connected (FC) layer followed by a sigmoid activation. Equation (9) details the entire computational procedure. Then, the predicted value, $\widehat{y}$ is obtained by taking the $\mathrm{sigmoid}$ function to the output of an $\mathrm{FC}$ layer as follows:

$$\widehat{y}=\mathrm{sigmoid}\left(\mathrm{FC}({Z_2}_d\parallel {Z_2}_c)\right),$$

(9)

where the ${Z_2}_d$ and ${Z_2}_c$ are referred to as the drug’s and cell’s embedding from the $Z_2$, respectively, and ∥ describes concatenation. Then the output y and ground truth $\widehat{y}$ are fed into the binary cross entropy loss L as below:

$$L(y,\widehat{y})=-{\displaystyle \frac{1}{N}}\sum _{i=1}^N\left[y_{i}log\left({\widehat{y}}_i\right)+\left(1-y_i\right)log\left(1-{\widehat{y}}_i\right)\right],$$

(10)

where N is the number of data. Note that this network utilized Adam [48] for optimization.

3.3. Model Configuration and Training

Our drGAT model was implemented using PyTorch 0.13.1 [49] and PyTorch Geometric 2.2.0 [50]. We utilized Optuna [51], a library that employs Bayesian Optimization, for hyperparameter tuning. Detailed hyperparameter settings can be found in Appendix B.5.

4. Discussion

In this study, we successfully developed the drGAT model, utilizing Graph Attention Networks to derive attention coefficients for each node, including genes, drugs, and cell lines. This approach is capable of binary drug response prediction, as well as, facilitating model interpretation. Our model demonstrated superior performance in terms of accuracy, recall, and F1 score, outperforming competing models in drug response prediction. Additionally, the drGAT model is capable of predicting drug sensitivity in untested drug-cell line combinations based on attention coefficients, achieving over 80% accuracy. This capability is valuable for suggesting drug efficacy to untested samples (evaluated with untested samples in from the training data).

A key highlight of our research is the ability of our drGAT model to predict drug-gene relationships beyond those included in the training data as drug-target interactions (DTI). For drug-gene relationships that are not DTIs, our analysis found many relationships that have been previously studied using the attention coefficient-based approach. This capability allows our model to not only produce accurate response predictions that utilize known drug-target interaction information but also identify potential drug-gene relationships. These relationships could be of further interest in understanding drug mechanisms. Additionally, as part of this work, we analyzed the attention coefficients to understand the biological processes relevant to each drug. Our analysis returned several key processes related to cancer, as well as, processes with genes relevant to cancer and drug response.

One challenge in this study was the modest volume of high-confidence drug-target interaction data available. Although the NCI60 dataset includes 269 DNA damage-related compounds, only 100 of the drug compounds had DrugBank drug-target interaction (for these drugs we included 571 interactions from DrugBank). This highlights the need for further collection of drug-gene relationships. Ongoing efforts by one of the present study’s authors and other colleagues are aimed at expanding such datasets. This expansion is actively utilizing advancements in information extraction from natural language processing as well as crowd-sourced curation efforts [6,52,53,54,55].

Overall, this work highlights how results from advanced machine-learning techniques can be paired with bioinformatic analyses to understand mechanisms of cellular response to treatment. We expect future work to broaden the approach while providing more granularity to its interpretability. Currently, our model leverages drug response data to generate cell line similarity. Future work could incorporate additional multi-omics data, including gene expression, methylation, copy number variation, and mutation data as has been done by others [56]. We believe this enhancement will boost both model performance and interpretability.

Code Availability: The implementation of the drGAT model, including the scripts used for training, validation, and evaluation, is available on GitHub https://github.com/inoue0426/drGAT. This repository contains all necessary dependencies, installation instructions, and guidelines for running the model.

Appendix B.1. Validation

The drug response matrix $X_{DR}$ is utilized as the input to the prediction model, using a split of 60% for training, 20% for validation, and 20% for testing. The drug response matrix contains z-score values that are transformed into a binary association matrix as follows: if the value in the drug response data for the NCI60 is positive, we consider that the cell line is sensitive to the drug and set the value to 1; otherwise, the cell line is resistant and set the value to 0.

Appendix B.2. Baseline Methods

For comparison, we have selected four recent studies as benchmark models. As discussed in subsection A, our work is related to MOFGCN, a current state-of-the-art method that uses a heterogeneous graph from drug response and gene expression. Consequently, MOFGCN is our primary benchmark model for GNN-based methods. Furthermore, we include a deep learning method to predict the drug sensitivity of cancer cell lines, called DeepDSC [73] for the DNN-based methods. DeepDSC employs a pipeline comprising six layers of Auto Encoder to reconstruct gene expression data. The hidden layer is extracted and subjected to convolution with the parameters of the Morgan fingerprint of the compounds. This convolved information is input for a DNN with four linear layers, predicting drug responses. The parameter settings for MOFGCN and DeepDSC were utilized from the original papers. We also include a DNN-based AutoML utilizing AutoKeras [74]. In addition, we utilize tree-based methods: LightGBM and Random Forest, tuning by FLAML [75].

Appendix B.3. Comparative GNN Architecture

Our study evaluates different Graph Neural Network (GNN) layers to determine their effectiveness. The layers compared are Message Passing Neural Networks (MPNN) [76], Graph Convolutional Networks (GCN) [71], Graph Attention Network (GAT) [17], Graph Attention Network V2 (GATv2) [46], and Graph Transformer (GT) [77]. GAT introduces an attention mechanism into graph neural networks. It allows nodes to assign varying levels of importance to their neighbors dynamically. GATv2 is an advancement over GAT. It features a more expressive attention mechanism that recalculates attention coefficients in every forward pass of the model. The GT layer extends the Transformer [16] architecture to graph-structured data. It leverages self-attention mechanisms specifically tailored for graphs.

Appendix B.4. Evaluation Metrics of Prediction Performance

To evaluate the prediction performance of our drGAT model, we employed four criteria: accuracy, precision, recall, and F1 scores. In binary classification, there are four kinds of test data points based on their ground truth and the model’s prediction,

• positive sample and is correctly predicted as positive, also known as True Positive (TP);
• negative samples and is wrongly predicted as positive samples, also known as False Positive (FP);
• negative samples and is correctly predicted as negative samples, also known as True Negative (TN);
• positive samples and is wrongly predicted as negative samples, also known as False Negative (FN).

These metrics are defined as

$$\begin{array}{cc}\hfill \mathrm{Accuracy}& ={\displaystyle \frac{\mathrm{TP}+\mathrm{TN}}{\mathrm{TP}+\mathrm{TN}+\mathrm{FP}+\mathrm{FN}}},\hfill \\ \hfill \mathrm{Recall}& ={\displaystyle \frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}},\hfill \\ \hfill \mathrm{Precision}& ={\displaystyle \frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FP}}},\hfill \\ \hfill \mathrm{F}1\mathrm{Score}& ={\displaystyle \frac{2\times \mathrm{Precision}\times \mathrm{Recall}}{\mathrm{Precision}+\mathrm{Recall}}},\hfill \end{array}$$

(A1)

where $\mathrm{TP}$, $\mathrm{FP}$, $\mathrm{TN}$, and $\mathrm{FN}$ correspondingly denote True Positive, False Positive, True Negative, and False Negative counts, respectively.

Appendix B.5. Hyperparameter Tuning

We utilized Optuna [51], a library that employs Bayesian Optimization to tune hyperparameters. Specifically, the number of epochs is 1500, the number of attention heads is 5, the number of hidden units is configured with 256 units for the first linear layer, 32 units for the second GNN layer, 128 units for the final GNN layer, the dropout ratio is 0.1, and the learning rate is 0.001. The total of parameters is 1,074,401.

Our experiments were conducted on an Ubuntu server equipped with an NVIDIA A40 GPU with 48 GB of memory. This hardware configuration allowed for efficient training and evaluation of our model, particularly beneficial for processing the large-scale heterogeneous graph data used in our study.