Graphical abstract graphic file with name fx1.jpg [40]Open in a new tab Highlights * • M3NetFlow is a graph AI model for integrative and interpretable multi-omic analysis * • M3NetFlow supports target and pathway inference based on given targets of interest * • M3NetFlow supports generic target and pathway inference from multi-omic data * • NetFlowVis can visualize multi-omic features of predicted targets and pathways __________________________________________________________________ Biocomputational method; Complex systems; Omics Introduction Multi-omic data-driven studies are at the forefront of precision medicine and healthcare. Recently, multi-omic datasets, like genetic, epigenetic, transcriptomic, and proteomic, have been generated to characterize dysfunctional biological processes and signaling pathways from multiple levels/views and to elucidate the panoramic view of the disease pathogenesis.[41]^1^,[42]^2^,[43]^3 For example, The Cancer Genome Atlas (TCGA) program has generated multi-omic datasets of over 20,000 samples spanning 33 cancer types, to understand the key molecular targets and signaling pathways of cancer. Moreover, the multi-omic data of >10,000 cancer cell lines were profiled in the Cancer Cell Line Encyclopedia (CCLE) project, which are valuable to investigate the mechanism of cancer response to given drugs and drug combinations.[44]^4 In addition, the multi-omic data of Alzheimer’s disease (AD) are generated and publicly available in The Religious Orders Study and Memory and Aging Project (ROSMAP[45]^5) project to uncover the pathogenesis of AD. Also, the exceptional longevity (EL), like the Long-Life Family Study project, has been generating multi-omic data[46]^6^,[47]^7 to identify protective biomarkers and pathways for long and healthy life. The multi-omic data are valuable and essential for understanding the key molecular targets and mechanisms of diseases, identifying novel therapeutic targets, predicting effective drugs and drug cocktails to guide the development of precision medicine. However, it remains an open problem and a challenging task to integrate multi-omic data and mine core disease signaling pathways from the complex and intensive signaling interactions among a large number of proteins in the cell signaling system.[48]^8 The two problems to be tackled in this study are (1) hypothesis/anchor-target guided multi-omic data analysis and (2) generic multi-omic data analysis. The hypothesis/anchor-targets can be known as disease-associated targets, or drug targets; and the expected outcome is the upstream or downstream signaling pathways of given anchor-targets. Specifically, the drug combination synergy score prediction task (as the hypothesis/anchor-target guided analysis example) is defined as follows. In two experimental datasets (National Cancer Institute [NCI] 60 and the O’Neil), the synergy scores (label information to train the model) of a set of pairwise drug combinations were experimentally measured on a set of cancer cell lines. The multi-omic data of these cancer cell lines were also measured and publicly accessible. The biological and computational problems are if a computational model can (1) predict synergy score of drug combinations and (2) investigate the potential mechanism of synergy (MoS) using (model inputs) the multi-omic data of cancer cell lines, known drug targets (drug-target interactions), and Kyoto Encyclopedia of Genes and Genomes (KEGG[49]^9) signaling pathways (a graph with signaling/protein-protein interactions). So each data point for the model can be represented as (in addition to KEGG signaling pathways). The expected model outcomes/output are (1) synergy score prediction capability of the model and (2) MoS of effective drug combination on given cell lines, i.e., important signaling targets and cascades/flows within cell lines that can explain which drug combinations are effective and which are not. For the generic multi-omic data analysis demonstrated using AD sample classification task, the multi-omic datasets of AD and control human samples are publicly accessible. The biological and computational problems are if a computational model can (1) predict sample categories and (2) identify disease-associated targets and pathways. Specifically, the model inputs are the (1) multi-omic data of samples, (2) KEGG signaling graph, and (3) sample categories (label information to train the model): AD vs. control. So each data point for the model can be represented as (in addition to KEGG signaling pathways). The expected model outputs are (1) model prediction capability to classify samples into AD vs. control and (2) the potential AD-associated targets. A comprehensive review of existing multi-omic data integration analysis models was reported.[50]^8 Specifically, these models were clustered into a few categories, like similarity, correlation, Bayesian, multivariate, fusion, and network-based models.[51]^8 The pathway representation and analysis by direct inference on graphical models (PARADIGM[52]^10) is one of the most widely used methods among these traditional computational methods. The mixOmics[53]^11 and timeOmics[54]^12 models and tools were developed to identify disease-associated genes using multivariate analysis on the multi-omic datasets of one and multiple time points, respectively, in which the signaling network information was not integrated. The mechanism of action generator involving network analysis (MAGINE)[55]^13 is an enrichment-based framework, in which the differentially expressed genes from multi-omic data are identified first, and enriched pathways and network modules are then identified based on the identified genes. The GLUE[56]^14 (graph-linked unified embedding) model was proposed to integrate multi-omic data by mapping raw omic data/features into a new embedding space for the cell clustering and correlation analysis. In the GLUE model, the signaling network information was converted to embedding features to adjust the embedding of raw omic data, which did not directly use the network information to calculate the signaling flow on the network. Recently, graph neural networks (GNNs) have gained prominence due to their capability to model relationships within graph-structured data.[57]^15^,[58]^16^,[59]^17 And numerous studies have applied the GNN with the integration of the multi-omic data. Multi-omics graph convolutional networks (MOGONET[60]^18) initially creates similarity graphs among samples by leveraging each omic data and then employs a graph convolutional network (GCN[61]^17) to learn a label distribution from each omic data independently. Subsequently, a cross-omic discovery tensor is implemented to refine the prediction by learning the dependency among multi-omic data. Multi-omics integration model based on graph convolutional network (MoGCN[62]^19) adopts a similar approach by constructing a patient similarity network using multi-omic data and then using GCN to predict the cancer subtype of patients. The universal framework for the integration of single-cell multi-omics data based on graph convolutional network (GCN-SC[63]^20) utilizes a GCN to combine single-cell multi-omic data derived from varying sequencing methodologies. Nevertheless, none of these models contemplate incorporating structured signaling data like KEGG into the model. Moreover, general GNN models are limited by their expression power, i.e., the low-pass filtering or over-smoothing issues, which hamper their ability to incorporate many layers. The over-smoothing problem was firstly mentioned by extending the propagation layers in GCN.[64]^21 Moreover, theoretical papers using Dirichlet energy showed diminished discriminative power by increasing the propagation layers. And multiple attempts were made to compare the expressive power of the GCNs,[65]^22 and it is shown that Weisfeiler-Lehman (WL) subtree kernel[66]^23 is insufficient for capturing the graph structure. Hence, to improve the expression powerful of GNN, the [MATH: K :MATH] -hop information of local substructure was considered in various recent research.[67]^24^,[68]^25^,[69]^26^,[70]^27^,[71]^28 However, none of these studies was specifically designed to well integrate the biological regulatory network and provide the interpretation with important edges and nodes. In this study, we present a novel graph AI model, M3NetFlow (multi-scale, multi-hop, multi-omic network flow) to address the challenges mentioned earlier. The major and unique contributions are as follows. Compared with existing models, we first mapped multi-omic data and drug-target information onto KEGG signaling pathways,[72]^9 which can support both anchor-target (drug targets) guided analysis and generic multi-omic analysis. To improve the model expression power, then we conducted the message propagation or signaling flowing with attention, multi-hop, on both local signaling modules and then updated the node embedding on the global signaling graph. This model design and model architecture are novel for multi-omic analysis tasks. Then the important signaling targets and interactions can be identified using the attention-based scores, learned in the model on the training dataset. To assess and demonstrate the effectiveness of the proposed model, M3NetFlow, it was applied in two independent multi-omic case studies: (1) uncovering mechanisms of synergy of effective drug combinations (hypothesis/anchor-target guided multi-omic analysis) and (2) identifying biomarkers of AD (generic multi-omic analysis). The evaluation and comparison results showed that M3NetFlow achieved the best prediction accuracy and identified a set of essential drug combination synergy- and disease-associated targets. To facilitate investigating the analysis results, a visualization tool, NetFlowVis, was developed to visualize the top-ranked signaling targets and interactions based on the attention-based target and interaction scores. The details of the studies are introduced in the following sections. Results [73]Figure 1 shows the schematic architecture of the proposed M3NetFlow model. The model input parameters are [MATH: X={(X(1< mo>),T(1< /mn>)),(X(2< mo>),T(2< /mn>)),…,(X(m< mo>),T(m< /mi>)),…,(X(M< mo>),T(M< /mi>))}(X(m< mo>)∈Rn×d< /mrow>,T∈Rn×2< /mrow>),A∈Rn×n< /mrow>,S={S1, S2,…,Sp,…,S P}(Sp∈Rnp×np),Din∈Rn×n< /mrow>,Dout< /mi>∈Rn×n< /mrow> :MATH] , where [MATH: M :MATH] represents the number of data points, [MATH: X :MATH] denotes all of the data points in the dataset, and [MATH: (X(m< mo>),T(m< /mi>)) :MATH] is [MATH: m :MATH] -th data points in the dataset, where [MATH: X(m) :MATH] denotes the node features matrix with [MATH: n :MATH] nodes of [MATH: d :MATH] features and [MATH: T(m) :MATH] denotes the one-hot encoding of two drugs (drug combinations) targeted on those [MATH: n :MATH] nodes (in the general multi-omic data analysis, the variable [MATH: T :MATH] will not be used). The matrix [MATH: A :MATH] is the adjacency matrix that demonstrates the node-node interactions. The element in adjacency matrix [MATH: A :MATH] such as [MATH: aij :MATH] indicates an edge from [MATH: i :MATH] to [MATH: j :MATH] . [MATH: S :MATH] is a set of subgraphs that partition the whole graph adjacent matrix [MATH: A :MATH] into multiple subgraphs with [MATH: Sp∈Rnp×np :MATH] of node interactions between its internal [MATH: np :MATH] nodes, and each subgraph has its own corresponding subgraph node feature matrix [MATH: Xp∈Rnp×d :MATH] . [MATH: Din :MATH] is an in-degree diagonal matrix for nodes in directed graph, and [MATH: Dout :MATH] is an out-degree diagonal matrix for nodes in directed graph. The model is to build up a model [MATH: f(·) :MATH] to predict the labels of samples: [MATH: f(X,A,S ,Din,Dout)=Y :MATH] , where [MATH: Y :MATH] is the sample labels, [MATH: Y∈RM×1< /mrow> :MATH] (e.g., the drug combination synergy scores or AD vs. control). Figure 1. [74]Figure 1 [75]Open in a new tab Model architecture of M3NetFlow (A) Mapping multi-omic data and drug targets onto KEGG signaling pathways. (B) Multi-hop attention-based signaling propagation on subgraphs. (C) Global signaling propagation. (D) Downstream tasks: (D.1) hypothesis/anchor-target guided decoder/prediction; (D.2) generic pooling/sorting based decoder/prediction. Experimental setup For drug combination synergy score prediction task, the 5-fold cross-validation ( [MATH: F=5 :MATH] ) was employed. For each sample, it has four elements: representing that drugs D[A] and D[B] are used on cell line C[C] and S[ABC] is the drug synergy score. The drug targets of D[A]and D[B] and the multi-omic data of cell line C[C] were used as the model input to predict S[ABC] (see [76]STAR Methods). Specifically, there are 2,788 samples for the NCI ALMANAC[77]^29 (A Large Matrix of Anti-Neoplastic Agent Combinations) drug combination dataset and 1,008 samples for the O’Neil dataset.[78]^30 For this task, 1,489 proteins on KEGG signaling graph were selected (overlapping with the omic data), and each protein (graph node) is characterized by 6 omic features: gene expression (RNA sequencing), copy-number variation (CNV), gene amplification, gene deletion, and gene methylation maximum and minimum values, as well as it is a target of D[A] or D[B] from Cell Model Passports,[79]^31 CCLE, and DrugBank[80]^32 databases. For AD sample classification task, multi-omic data of 138 (74 AD, 64 control [non-AD]) samples were derived from the ROSMAP[81]^5 database. We randomly selected 64 AD and 64 control samples as a balanced dataset and used the 5-fold cross-validation ( [MATH: F=5 :MATH] ) to evaluate the model performance. For this task, 2,099 proteins on KEGG signaling graph were selected (overlapping with the omic data), and each protein (graph node) is characterized by 10 omic features: methylation values on upstream, distal promoters, proximal promoters, core promoters and downstream, genetic duplication, deletions, CNV, and gene expression and protein expression from ROSMAP and GEO[82]^33 platform. Hyperparameters The models were implemented using pytorch and torch geometric. For both tasks, the learning rate started at 0.002 and was reduced equally within each batch for a certain epoch stage. After 60 epochs, the learning rate was set at 0.0001. Adam optimizer was used for optimization with eps = 1e−7 and weight_decay = 1e−20. We empirically set the [MATH: K :MATH] -hop subgraph message propagation with [MATH: K=3 :MATH] (3 hops) and the global bi-directional message propagation with [MATH: L=3(3layers) :MATH] . Afterward, the feature dimensions will vary at the different layers and will be denoted by [MATH: (d(1< mo>),d(2< /mn>),…,< mi>d(l), …,d(L )) :MATH] . At the global message propagation step, the layer will concatenate new node embedding features and output node features of the previous layer for both upstream and downstream, generating concatenated dimensions for the output dims being [MATH: 3×d(l ) :MATH] in [MATH: l :MATH] -th layer. Output dims of the previous layer served as the input dims for the current layer, as follows: (1) first layer’s (input dims, output dims), [MATH: (d(0< mo>),3d(< /mo>1)) :MATH] ; (2) second layer’s (input dims, output dims), [MATH: (3d(< mn>1),3d< /mi>(2)) :MATH] ; and (3) third layer’s (input dims, output dims), [MATH: (3d(< mn>2),3d< /mi>(3)) :MATH] . The final embedded node dims were [MATH: 3d(3 )(L=3) :MATH] . For drug combination synergy score prediction task, the decoder trainable transformation matrix dims [MATH: D∈R3d(L)×E :MATH] and [MATH: U∈RE×E< /mrow> :MATH] were used as trainable decoder matrices; [MATH: E=150 :MATH] was used. In AD sample classification task, the max pooling was employed to predict the sample outcome as described in [83]Equation 6. As for the LeakyReLU function, the parameter was set as 0.1 for both tasks. M3NetFlow improves prediction accuracy To evaluate the model performance in terms of synergy score prediction for drug combinations and predictions on ROSMAP AD samples, we conducted 5-fold cross-validation. As shown in [84]Table 1, the average prediction (using the Pearson correlation coefficient) was about 61% Pearson correlation using the test data in the NCI ALMANAC dataset and was about 64% Pearson correlation using the test data in the O’Neil dataset. Regarding the ROSMAP dataset, the average prediction accuracy was about 66% using the test data in the ROSMAP dataset. These prediction results are comparable with existing deep learning models.[85]^34^,[86]^35 Moreover, we also compared our proposed model M3NetFlow with other deep learning models, which included the GCN,[87]^17 graph attention network[88]^16 (GAT), UniMP,[89]^36 MixHop,[90]^25 principal neighborhood aggregation[91]^37 (PNA), and graph isomorphism network (GIN). By checking the p values over 5-fold cross-validation, the performances of the M3NetFlow have significant improvement over most of the GNN-based methods (see [92]Table 1 and [93]Figure 2A). Table 1. Model comparisons using average Pearson correlation and prediction accuracy (mean ± standard deviation) of 5-fold cross-validation using NCI ALMANAC, O’Neil, and ROSMAP datasets Dataset NCI ALMANAC O’Neil ROSMAP GCN 51.93% ± 3.55% 44.47% ± 6.60% 59.43% ± 4.53% GAT 49.16% ± 2.26% 57.06% ± 3.72% 62.80% ± 7.11% UniMP 49.02% ± 4.62% 55.84% ± 10.93% 61.83% ± 3.78% MixHop 57.78% ± 3.66% 27.15% ± 10.01% 57.20% ± 3.92% PNA 55.63% ± 2.47% 62.20% ± 2.27% 57.83% ± 1.82% GIN 53.76% ± 2.47% 33.12% ± 9.89% 49.83% ± 5.71% M3NetFlow 60.72% ± 0.77% 64.36% ± 2.53% 67.34% ± 5.12% [94]Open in a new tab Figure 2. [95]Figure 2 [96]Open in a new tab Model performance and overview of input datasets NCI ALMANAC, O’Neil (drug combination multi-omic data), and ROSMAP (AD multi-omic data) (A) Averaged Pearson correlation across 5-fold validation comparisons for GCN, GAT, UniMP, MixHop, PNA, GIN, and M3NetFlow models for NCI ALMANAC, O’Neil, and ROSMAP dataset (data are represented as mean). (B) Scatterplot of the model with data points in the whole NCI ALMANAC dataset. (C) Distributions of all cell lines in the whole NCI ALMANAC dataset. (D) Boxplots across all cell lines in the whole NCI ALMANAC dataset. M3NetFlow ranks important targets and interactions via attention scores Targets with higher importance scores predicting drug combination synergy Based on the node importance score calculated from edge attention scores (see [97]target/node importance score calculation in [98]STAR Methods section for details), the important targets for each cell line can be selected. For example, to investigate the MoS, we compared the importance scores of drug targets of the top 5 drug combinations (with highest synergy scores) and bottom 5 drug combinations (lowest drug synergy scores) ([99]Figure 3A). In another word, it is expected that the targets of synergistic drug combinations have higher importance scores than the targets of the non-synergistic drug combinations. Our results confirmed this pattern in 37 out of 41 (∼90%) cell lines, which have higher mean target importance scores (light green in [100]Table S1). Specifically, in 27 out of 41 (∼65%) cell lines, the targets of synergistic drug combinations have higher importance scores with p value ≤ 0.1 (deep orange color in [101]Table S1), and, in 32 out of 41 (∼78%) cell lines, the targets of synergistic drug combinations have higher importance scores with p value ≤ 0.3 (light orange color in [102]Table S1). [103]Figure 4 shows the pattern of target/node importance scores of drug combinations. Specifically, the left boxplots compared the node/target importance score distribution of individual top 5 and bottom (low) 5 drug combinations, respectively (each drug combination has multiple targets). The right boxplots show the node/target importance score distribution of all top 5 and bottom 5 drug combinations, respectively. The scores of individual drug targets and proteins in each cell line were provided in [104]Table S2, and the top 20 protein lists for each cell line were provided in [105]Table S3. Moreover, we identified the commonly important proteins across cell lines of the same cancer type ([106]Figure S2). Figure 3. [107]Figure 3 [108]Open in a new tab Target importance score patterns of synergistic and non-synergistic drug combinations (A) Illustration of analyzing the important targets of top 5 synergistic and bottom 5 non-synergistic drug combinations. (B) Visualization tool NetFlowVis for core signaling network interactions of cell line DU-145. (C–F) Target importance score distribution and boxplots of the top 5 synergistic and bottom 5 non-synergistic drug combinations in DU-145 and SK-MEL-28 cell lines, and t-test p-values were used to statistically compare the two distributions in (C) and (E). Figure 4. [109]Figure 4 [110]Open in a new tab Boxplots of target importance scores of cell lines A498, A549/ATCC, ACHN, BT-549, CAKI-1, DU-145, EKVX, HCT-116, HCT-15, HOP-62, HOP-92, HS 578T, IGROV1, K-562, KM12, LOX IMVI, MCF7, MDA-MB-231/ATCC, MDA-MB-468, NCI-H23, HCI-H460, NCI-H522, OVCAR-3, OVCAR-4, OVCAR-8, PC-3, RPMI-8226, SF-268, SF-295, SF-539, SK-MEL-28, SK-MEL-5, SK-OV-3, SNB-75, SR, SW-620, T-47D, U251, UACC-257, UACC-62, and UO-31 Top-ranked targets are associated with AD By setting the filters (edge threshold as 0.106 and a small component threshold as 15), we identified 100 potential important genes for AD. Among those genes, 28 genes are filtered out by setting the threshold of attention-based node weight as 2.0 and 15 of them with p values smaller than 0.1 in at least one of the 10 multi-omic features (see [111]Figures 5A and 5B). To evaluate the top targets ranked by attention-based node weight, the top-ranked 28 AD-associated biomarkers were further analyzed via pathway enrichment analysis. Interestingly, a set of AD-associated signaling pathways are identified. As shown in [112]Figure 5C, the top-ranked targets are involved in a set of signaling transduction pathways, which indicates the importance of these targets. Among them, several signaling pathways play critical roles in AD particularly through mechanisms involving inflammation, immune responses, and cellular growth and death. For instances, the B cell receptor (BCR) and T cell receptor (TCR) signaling pathways are vital for B cell and T cell activation, which plays a crucial role in immune surveillance and inflammation. Shared molecular components between the BCR and TCR pathways, including Mitogen-Activated Protein Kinase Kinase 1 (MAP2K1), Jun Proto-Oncogene, AP-1 Transcription Factor Subunit (JUN), Component Of Inhibitor Of Nuclear Factor Kappa B Kinase Complex (CHUK), RELA Proto-Oncogene, NF-KB Subunit (RELA), Nuclear Factor Kappa B Subunit 1 (NFKB1), Inhibitor Of Nuclear Factor Kappa B Kinase Subunit Beta (IKBKB), and protein kinase B (AKT) genes, suggest significant crosstalk that contributes to neuroinflammatory processes in AD. Prolonged activation of these immune pathways may exacerbate chronic inflammation and contribute to AD pathology by sustaining harmful neuroinflammatory responses.[113]^38 The nuclear factor κBNF-κB signaling axis, a critical component in both BCR and TCR pathways, is particularly relevant in AD due to its role in regulating inflammatory cytokine production. Chronic activation of NF-κB has been linked to increased amyloid-beta (Aβ) deposition and neurofibrillary tangle formation, both of which are hallmark features of AD pathology.[114]^39^,[115]^40 Additionally, activation of NF-κB, mitogen-activated protein kinase (MAPK), and AKT signaling cascades can induce neuronal apoptosis, leading to cognitive decline associated with AD. The MAPK, rat sarcoma (RAS), Phosphatidylinositol 3-Kinase / Protein Kinase B (PI3K/Akt), and forkhead box O (FoxO) signaling pathways also play critical roles in modulating immune responses and neuronal survival. The p38 MAPK pathway, for example, drives neuroinflammation by activating microglia and promoting the release of pro-inflammatory cytokines, which can result in neuronal death.[116]^41 Overactivation of RAS signaling increases oxidative stress, contributing to Aβ accumulation and tau pathology, both of which are central to neurodegeneration in AD.[117]^42 While Aβ accumulation is an early event in AD, tau pathology is more closely associated with cognitive decline, and, together, these processes are considered the primary drivers of neuronal death in AD.[118]^43^,[119]^44 The Akt signaling pathway, which promotes cell survival by inhibiting apoptosis through downstream effectors such as B-cell CLL/lymphoma 2 (BCL2), is also impaired in AD. Reduced Akt activity has been linked to increased neuronal apoptosis and the accumulation of Aβ and hyperphosphorylated tau. Notably, the PI3K-Akt pathway is intimately connected to insulin signaling, which is disrupted in AD and is characterized by reduced Akt signaling, impaired glucose metabolism, and an increased vulnerability to neurodegeneration.[120]^45 Furthermore, neurotrophin signaling, which regulates axonal growth and regeneration via MAPK and PI3K-Akt pathways, is another pathway that becomes disrupted in AD. This disruption impairs axonal repair mechanisms and contributes to synaptic loss and neuronal degeneration.[121]^46 Additionally, endocrine-related pathways, including insulin and relaxin signaling, play crucial roles in maintaining cellular homeostasis and survival. Dysregulation of these pathways in AD contributes to neuronal apoptosis, neuroinflammation, and impaired protein clearance, further exacerbating disease pathology.[122]^45^,[123]^47^,[124]^48 Figure 5. [125]Figure 5 [126]Open in a new tab Top-ranked proteins and associated pathways related to Alzheimer’s disease (A) Visualization tool NetFlowVis-AD for core signaling network interactions of AD. Red color indicates nodes with node score > 2.0. The purple circle represents nodes with at least one feature’s p value < 0.1 between AD and control samples. (B) Barplot of the node scores, where the red bars and purple borders have the same meaning as in (A). The color indicates the p value of individual omic data between AD and control samples. (C) Sankey plot of top-enriched signaling pathways of the top-selected proteins. Discussion Multi-omic data-driven studies are the forefront of biomedical research. Large-scale multi-omic datasets have been generated to characterize the dysfunctional targets and signaling pathways of complex diseases, like cancer and AD, which are valuable and essential for the development of personalized medicine or precision medicine. However, it remains an open problem to integrate and interpret the multi-omic datasets to identify key molecular targets and core signaling pathways. In this study, we clearly formulated and defined the two types of need of multi-omic data analysis, i.e., hypothesis/anchor-target guided multi-omic data analysis and generic multi-omic data analysis. To tackle the multi-omic data analysis tasks, we proposed a graph AI model, M3NetFlow, which is specifically designed for integrating and analyzing multi-omic datasets and can deal with both of analysis tasks. Compared with existing models, the proposed model design and model architecture, mapping multi-omic data into signaling networks, combing local and global signaling, and multi-hop flowing/propagation on the signaling network, are unique. The model prediction and interpretation (attention-based target and interaction scores) capabilities are demonstrated and evaluated using two case studies: drug combination synergy score prediction task (hypothesis/anchor-target guided multi-omic data analysis) and AD sample classification (generic multi-omic data analysis). The source code of M3NetFlow is publicly accessible, which enables users to modify and improve the model for their own studies. This method can be an alternative option and can be combined with other analysis methods, like the MAGINE for pathway and network module enrichment analysis. Limitations of the study This study is an exploratory effort in multi-omic data analysis, with several areas requiring further investigation. For example, more signaling pathways and larger protein-protein interactions should be evaluated. Moreover, dividing large signaling graphs into subnetworks or network modules can be achieved by using biologically meaningful annotations, such as gene ontology (GO) terms. In the current analysis, the generic pre-analyzed multi-omic features were used. It is worth investigating and selecting additional and biological meaningful features that can be derived from the multi-omic datasets in the future work. Additionally, more and more multi-omic datasets are being generated. Combining multi-omic data from different diseases can provide a larger sample size than individual disease datasets, which could improve the training or pre-training of graph AI models and help identify pan-disease or disease-specific targets. It is also interesting to expand graph models from tissue-level multi-omic data to single-cell multi-omic data, which can be more challenging due to the large number of single-cell samples. Aside from this, incorporating large language models (LLMs) into graph-based frameworks offers further possibilities, such as interpreting textual annotations/descriptions and associated GO terms or pathways to enrich the informative features, which can improve prediction accuracy and identify the accurate and interpretable biomarkers and core signaling pathways. Developing a hybrid graph-language framework that integrates structured data with unstructured text could enable contextual learning, while pre-trained LLMs like Generative Pre-training Transformer (GPT) or Bidirectional Encoder Representations from Transformer (BERT) may assist in mapping textual annotations to graph components, improving pathway discovery. Therefore, sophisticated and improved graph AI models are needed to integrate and interpret multi-omic datasets, identify and infer key molecular targets and signaling pathways of complex diseases, and guide the development of precision medicine. Resource availability Lead contact Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact Prof. Fuhai Li (fuhai.li@wustl.edu). Materials availability This study did not generate unique reagents. Data and code availability * • The multi-omic data derived from cancer cell lines and AD samples, along with the drug synergy dataset and the dataset detailing AD sample conditions, are accessible through the [127]key resources table. * • The data processing code and associated details, along with the M3NetFlow code, are available in the GitHub repository links given in the [128]key resources table. Acknowledgments