Pathformer: a biological pathway informed transformer for disease diagnosis and prognosis using multi-omics data

Authors: Xiaofan Liu, Yuhuan Tao, Zilin Cai, Pengfei Bao ⚬, Hongli Ma, Kexing Li, Mengtao Li ⚬, Yunping Zhu, Zhi John Lu ⚬

Affiliations: MOE Key Laboratory of Bioinformatics, Center for Synthetic and Systems Biology, School of Life Sciences, Tsinghua University, Beijing 100084, China; Institute for Precision Medicine, Tsinghua University, Beijing 100084, China; Department of Rheumatology and Clinical Immunology, Peking Union Medical College Hospital, Chinese Academy of Medical Sciences, Peking Union Medical College, National Clinical Research Center for Dermatologic and Immunologic Diseases (NCRC-DID), MST State Key Laboratory of Complex Severe and Rare Diseases, MOE Key Laboratory of Rheumatology and Clinical Immunology, Beijing 100730, China; State Key Laboratory of Medical Proteomics, Beijing Proteome Research Center, National Center for Protein Sciences (Beijing), Beijing Institute of Lifeomics, Beijing 102206, China

License: © The Author(s) 2024. Published by Oxford University Press. CC BY 4.0 This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.

Article links: DOI: 10.1093/bioinformatics/btae316  |  PubMed: 38741230  |  PMC: PMC11139513

Relevance: Moderate: mentioned 3+ times in text

Full text: PDF (7.1 MB)

1 Introduction

Comparing to a single type of data, multi-omics data provide a more comprehensive view of gene regulation (ref. Hasin ). Therefore, integrating multi-omics data from tissue and liquid biopsies would be helpful in addressing challenges in disease diagnosis (ref. Ning ), treatment (ref. Chiu , ref. Sharifi-Noghabi ), and prognosis (ref. Hao ), such as deregulated network between different types of molecules and data noise caused by patients’ heterogeneity (ref. Tarazona ). To integrate multi-omics data of cancer, several supervised methods have been developed, such as mixOmics (ref. Rohart ), liNN (ref. Kuru ), eiNN (ref. Preuer ), liCNN (ref. Islam ), eiCNN (ref. Fu ), MOGONet (ref. Wang ), and MOGAT (ref. Xing ). Later, the performance and interpretability of multi-omics data integration were further improved using deep learning models informed by biological pathways. For instance, a pathway-associated sparse deep neural network (PASNet) was utilized to predict the prognosis of glioblastoma multiforme (GBM) patients (ref. Hao ). Recently, P-NET, a sparse neural network integrating multiple molecular features based on a multilevel view of biological pathways, was introduced to predict subtype and survival of prostate cancer patients (ref. Elmarakeby ). In addition, PathCNN based on a convolutional neural network (CNN) was developed to predict the prognosis of GBM patients using principal component analysis (PCA) to define image-like multi-omics pathways (ref. Oh ).

These pathway-informed deep learning methods did not consider the crosstalk between omics and between pathways, although the crosstalk holds biological significance as well as pathway itself (ref. Kim , ref. Li , ref. Prahallad and Bernards 2016, ref. Liu ). The crosstalk means a member of one pathway regulates a component of another pathway. The balance and oscillation between different pathways contribute to cancer progression and metastasis. For instance, a positive feedback loop between Wnt pathway and ERK pathway were revealed in cancer(ref. Kim ); crosstalk between TGF-$\beta$ pathway and TNF-$\alpha$ pathway can promote tumor’s invasion and metastasis by affecting its microenvironment (ref. Liu ). Meanwhile, the criss-cross attention mechanism of the Transformer would be very useful to capture the crosstalk information (ref. Jumper ). However, incorporating multi-omics data and their crosstalk information in a Transformer is very challenging: when processing multi-omics data, the multi-modal features are usually multiplied by tens of thousands of genes, producing an extremely long input that is not acceptable by a common Transformer model (usually <512 words). Meanwhile, certain embedding methods for biological data, such as discretization and linear transformation, were introduced in the previous Transformer models (ref. Osseni , ref. Cui , ref. Theodoris ), while biological information was largely lost during these kinds of embedding.

In order to integrate multi-omics data by embedding biological pathway crosstalk without information loss, we introduce a Transformer model, Pathformer, with three key steps to address the above problems. First, it transforms various modalities into distinct gene-level features using a series of statistical methods, such as the maximum value method, and connects these features into a novel compacted multi-modal vector for each gene, which not only preserves valuable information but also shortens the input. Second, Pathformer utilizes a sparse neural network based on prior pathway knowledge to transform gene embeddings into pathway embeddings. Third, Pathformer naturally incorporates pathway crosstalk network into a Transformer model with bias to enhance the exchange of information between different pathways and between different modalities (e.g. omics) as well.

Here, we first benchmarked Pathformer and 18 other integration methods in various classification tasks, using multiple cancer tissue datasets from TCGA. Then, we used Pathformer to integrate various multi-omics data from tissue and liquid biopsies. Through case studies on survival prediction of breast cancer and noninvasive diagnosis of pan-cancer, we revealed interesting pathways, genes, and regulatory mechanisms related to cancer in human tissue and plasma, demonstrating the prediction accuracy and biological interpretability of Pathformer in various clinical applications.

2 Materials and methods

2.1 Overview of Pathformer

Pathformer is mainly designed to integrate various multi-omics data from tissue and liquid biopsies, which can be used for different classification tasks in disease diagnosis and prognosis, such as cancer early detection, cancer staging and survival prediction (Fig. 1a). It has six modules: (i) biological pathway and crosstalk network calculation module, (ii) multi-omics data input module, (iii) biological multi-modal embedding module (key module), (iv) Transformer module with pathway crosstalk network bias, (v) classification module, and (vi) biological interpretability module.

Pathformer combines prior biological pathway information (module 1, Fig. 1a) with multi-modal data (module 2, Fig. 1b) for disease diagnosis and prognosis. It introduces a new embedding method to incorporate biological multi-modal data at both gene level and pathway level: it initiates the process by uniformly transforming different modalities to the gene level through a series of statistical indicators, then concatenates these modalities into compacted multi-modal vectors to define gene embedding, and employs a sparse neural network based on the gene-to-pathway mapping to transform gene embedding into pathway embedding (module 3, Fig. 1c). Pathformer then enhances the fusion of information between various modalities and pathways by combining pathway crosstalk networks with Transformer encoder (module 4, Fig. 1a, Supplementary Fig. S1). Finally, a fully connected layer serves as the classifier for different downstream classification tasks (module 5). In addition, Pathformer uses a biological interpretable module with attention weights and SHapley Additive exPlanations (ref. Lundberg and Lee 2017) values to identify important genes, pathways, modalities, and their crosstalk or regulation (module 6). These six modules are described in detail below.

2.2 Module 1: curation of biological pathways and calculation of initial crosstalk network

We curated 2289 biological pathways from four public databases including Kyoto Encyclopedia of Genes and Genomes database (KEGG) (ref. Kanehisa and Goto 2000), Pathway Interaction database (PID) (ref. Schaefer ), Reactome database (Reactome) (ref. Croft ), and BioCarta Pathways database (BioCarta) (ref. Nishimura 2001). Then, we filtered these pathways by three criteria: the gene number, the overlap ratio with other pathways (the proportion of genes in the pathway that are also present in other pathways), and the number of pathway subsets (the number of sub-pathways included in the pathway). Following the principle of moderate size and minimal overlap with other pathway information, we selected 1497 pathways with gene number between 15 and 100, or gene number >15 and overlap ratio <1, or gene number >15 and the number of pathway subsets <5. Next, we used BinoX (ref. Ogris ) to calculate the crosstalk relationship of pathways and build a pathway crosstalk network with adjacency matrix $\mathit{P}\in {\mathbb{R}}^{N_p{\times N}_p}$, $N_p$ = 1497 (more details in Supplementary Note S1).

2.3 Modules 2 and 3: multi-omics data input and multi-modal embedding

Biological multi-modal data preprocessing and embedding method are two key modules of Pathformer (Fig. 1b and c). In module 2 (Fig. 1b), to capture more comprehensive regulatory information, we expanded biological multi-omics data into multi-modal data, including not only data from different omics sources but also variant features of the same omics, such as RNA splicing, RNA editing, RNA alternative promoter, and so on. To obtain multi-modal data, we used standardized bioinformatics pipeline to calculate different omics or variant features of the same omics from raw sequence reads (more details in Supplementary Note S2). These multi-modal data have different dimensions, including nucleotide level, fragment level, and gene level. For example, Pathformer’s input for cancer tissue datasets from Cancer Genome Atlas (TCGA) (ref. Cancer Genome Atlas Research Network 2013) includes gene-level RNA expression, fragment-level DNA methylation, and both fragment-level and gene-level DNA CNV. Modalities and their dimension levels for different datasets are described in Supplementary Table S1.

In module 3 (Fig. 1c), we proposed a new biological multi-modal embedding method of Pathformer, which consists of gene embedding ${\mathit{E}}_G$ and pathway embedding ${\mathit{E}}_P$. We represented biological multi-modal input matrix of a sample as $\mathit{M}$, described as follows:

where m is the number of modalities, and $N_{M_i}$ is the length of input for modality i, like the number of genes for RNA expression, the number of editing sites for RNA editing, and the number of CpG islands for DNA methylation.

Next, we first used a series of statistical indicators to convert different modalities into gene level modal features, and then concatenated these modal features into a compressed multi-modal vector as gene embedding ${\mathit{E}}_G\in {\mathbb{R}}^{N_g\times D_g}$, which are calculated as follows:

where ${\mathit{G}}_l=[{\mathit{V}}_1^{\left(l\right)},\dots ,{\mathit{V}}_i^{\left(l\right)},\dots ,{\mathit{V}}_m^{\left(l\right)}]{\in \mathbb{R}}^{D_g}$ is gene embedding of the lth gene and is a compacted multi-modal vector; ${\mathit{V}}_i$ is the modal feature matrix of modality i; ${\mathit{V}}_i^j$ is the jth dimension of modal feature matrix for modality i; $N_g$ is the number of genes, $D_g=e_1+e_2+\cdots +e_m$ is the dimension of gene embedding; $e_i$ is the dimension of modal feature matrix for modality i; ${\mathit{F}}_E$ is the conversion function, which uses statistical indicators to uniformly convert different modalities into gene level; ${\mathit{F}}_{E_i}$ is the conversion function of modality i, and each modality’s function is constructed from distinct statistical indicator functions $f_{*}$ (more details in Supplementary Table S1). These statistical indicator functions include gene level score ($f_1$), count ($f_2$), minimum ($f_3$), maximum ($f_4$), mean ($f_5$), entropy ($f_6$), weighted mean in whole gene ($f_7$), and weighted mean in window ($f_8$), formulas of which are in Supplementary Note S3.

Subsequently, we used the known gene-pathway mapping relationship to develop a sparse neural network based on prior pathway knowledge (PSNN) to transform gene embedding ${\mathit{E}}_G$ into pathway embedding ${\mathit{E}}_P$, as described below:

where $N_p$ is the number of pathways, $D_p=D_g$ is the dimension of pathway embedding, ${\mathit{W}}_{\mathit{sparse}}\in {\mathbb{R}}^{N_g\times N_p}$ is a learnable sparse weight matrix, and $\mathit{B}$ is a bias term. ${\mathit{W}}_{\mathit{sparse}}$ is constructed based on the known relationship between pathways and genes. When the given gene and the pathway are irrelevant, the corresponding element of ${\mathit{W}}_{\mathit{sparse}}$ will always be 0. Otherwise, it needs to be learned through training. Therefore, pathway embedding is a dynamic embedding method. The PSNN cannot only restore the mapping relationship between genes and pathways, but also capture the different roles of different genes in pathways, and can preserve the complementarity of different modalities. Additionally, this biological multi-modal embedding step does not require additional gene selection, thereby avoiding bias and overfitting problems resulting from artificial feature selection.

2.4 Module 4: transformer module with pathway crosstalk network bias

We developed the Transformer module based on criss-cross attention (CC-attention) with bias for data fusion of pathways, modalities, and their crosstalk (Supplementary Fig. S1). This module has 3 blocks, each containing multi-head column-wise self-attention (col-attention), multi-head row-wise self-attention (row-attention), layer normalization, GELU activation, residual connection, and network update. Particularly, col-attention is used to enhance the exchange of information between pathways, with the pathway crosstalk network matrix serving as the bias for col-attention to guide the flow of information. Row-attention is employed to facilitate information exchange between different modalities, and the updated pathway embedding matrix is used to update the pathway crosstalk network matrix by calculating the correlation between pathways.

Multi-head column-wise self-attention contains 8 heads, and each head is a mapping of ${\mathit{Q}}_1$,${\mathit{K}}_1$,${\mathit{V}}_1$,$\mathit{P}$, which are query vector, key vector, and value vector of pathway embedding ${\mathit{E}}_P$ and pathway crosstalk network matrix $\mathit{P}$, respectively. First, we represented the hth column-wise self-attention by ${\mathit{A}}_{\mathit{col}}^{\left(h\right)}$, calculated as follows:

where $h=1,2,\dots ,H$ is the hth head; H is the number of heads; ${\mathit{W}}_{Q_1}^{\left(h\right)}\in {\mathbb{R}}^{D_p\times d}$, ${\mathit{W}}_{K_1}^{\left(h\right)}\in {\mathbb{R}}^{D_p\times d}$, ${\mathit{W}}_{V_1}^{\left(h\right)}\in {\mathbb{R}}^{D_p\times d}$ are the weight matrices as parameters; d is the attention dimension; ${\mathrm{dropout}}_{0.2}$ is a dropout neural network layer with a probability of 0.2; and softmax is the normalized exponential function.

Next, we merged multi-head column-wise self-attention and performed a series of operations as follows:

where $h=1,2,\dots ,H$ is the hth head; H is the number of heads; $\circ$ is the matrix dot product operator; ${\mathit{W}}_{g_1}^{\left(h\right)}\in {\mathbb{R}}^{D_p\times d}$, ${\mathit{W}}_{U_1}^{\left(h\right)}\in {\mathbb{R}}^{d\times D_p}$, ${\mathit{W}}_{O_{11}}\in {\mathbb{R}}^{D_p\times o}$, ${\mathit{W}}_{O_{12}}\in {\mathbb{R}}^{o\times D_p}$ are the weight matrices as parameters; o is a constant; LN is the layer normalization function; GELU is the distortion of RELU activation function; and ${\mathrm{dropout}}_{0.2}$is a dropout neural network layer with a probability of 0.2.

Multi-head row-wise self-attention enables information exchange between different modalities. It is a regular dot-product attention. It also contains eight heads, and the hth row-wise self-attention, i.e. ${\mathit{A}}_{\mathit{row}}^{\left(h\right)}$, is calculated as follows:

where $h=1,2,\dots ,$ h is the hth head; H is the number of heads; ${\mathit{W}}_{Q_2}^{\left(h\right)}\in {\mathbb{R}}^{N_p\times d}$, ${\mathit{W}}_{K_2}^{\left(h\right)}\in {\mathbb{R}}^{N_p\times d}$, ${\mathit{W}}_{V_2}^{\left(h\right)}\in {\mathbb{R}}^{N_p\times d}$ are the weight matrices as parameters; d is the attention dimension; ${\mathrm{dropout}}_{0.2}$ is a dropout neural network layer with a probability of 0.2; and softmax is the normalized exponential function.

Subsequently, we merged multi-head row-wise self-attention and performed a series of operations. The formulas are as follows:

where $h=1,2,\dots ,$ h is the hth head; H is the number of heads; $\circ$ is the matrix dot product operator; ${\mathit{W}}_{g_2}^{\left(h\right)}\in {\mathbb{R}}^{N_p\times d}$, ${\mathit{W}}_{U_2}^{\left(h\right)}\in {\mathbb{R}}^{d\times N_p}$, ${\mathit{W}}_{O_{21}}\in {\mathbb{R}}^{N_p\times o}$, ${\mathit{W}}_{O_{22}}\in {\mathbb{R}}^{o\times N_p}$ are the weight matrices as parameters; o is a constant; $\mathit{\beta }$ is a constant coefficient for row-attention; LayerNorm is the layer normalization function; GELU is the distortion of RELU activation function; and ${\mathit{dropout}}_{0.2}$ is a dropout neural network layer with a probability of 0.2. ${\mathit{O}}_2$ is pathway embedding input of the next Transformer block. In other words, when ${\mathit{E}}_P$ is ${\mathit{E}}_P^{\left(0\right)}$, ${\mathit{O}}_2$ is ${\mathit{E}}_P^{\left(1\right)}$. Superscripts with parenthesis represent data at different block.

Then, we used the updated pathway embedding ${\mathit{O}}_2$ to update the pathway crosstalk network. We exploited the correlation between embedding vectors of two pathways to update the corresponding element of the pathway crosstalk network matrix. The formula is as follows:

where ${\mathit{P}}^{\prime }$ is the updated pathway crosstalk network matrix of the next Transformer block. In other words, when ${\mathit{P}}^{\prime }$ is ${\mathit{P}}^{\left(1\right)}$, $\mathit{P}$ is ${\mathit{P}}^0$. Superscripts with parenthesis represent data at different block.

2.5 Module 5: classification module

Given the classification tasks in disease diagnosis and prognosis, we used the fully connected neural network as the classification module to transform pathway embedding encoded by the Transformer module into the probability for each label. Three fully connected neural networks each have 300, 200, and 100 neurons, with dropout probability ${\mathit{dropout}}_c$, which is a hyperparameter. More details are described in Supplementary Note S4.

2.6 Module 6: biological interpretability module

The biological interpretable module enables us to calculate the contribution of each modality, identify important pathways and their key genes, and uncover the most critical pathway crosstalk subnetworks.

To calculate the contribution of each omics and each modality, we first integrated all matrices of row-attention maps into one matrix by element-wise averaging. Then, we averaged this average row-attention matrix along with columns as the attention weights of modalities. More details are described in Supplementary Note S5.

To identify important pathways and their key genes, we used SHapley Additive exPlanations (ref. Lundberg and Lee 2017) (SHAP value) to calculate the contribution of each feature. It is an additive explanation model inspired by coalitional game theory, which regards all features as “contributors.” SHAP value is the value assigned to each feature, which explains the relationship between modalities, pathways, genes and classification, implemented by “SHAP” package of Python v3.6.9. Then, pathways with the top 15 SHAP values in the classification task are considered as important pathways. For each pathway, genes with top five SHAP values are considered as the key genes of the pathway. The modality of a gene with the rank of SHAP value higher than other modalities is considered the core modality of the gene. More details are described in Supplementary Note S5.

Particularly, the pathway crosstalk network matrix is used to guide the direction of information flow, and updated according to updated pathway embedding in each Transformer block. Therefore, the updated pathway crosstalk network contains not only the prior information in the initial network (module 1) but also the multi-modal data information derived from the Transformer module (module 4), which represents the specific regulatory mechanism in each classification task. We defined the sub-network score through SHAP value of each pathway in the sub-network, so as to find foremost sub-network for prediction, i.e. hub module of the updated pathway crosstalk network. The calculation of the sub-network score can be divided into four steps: average pathway crosstalk network matrix calculation, network pruning, sub-network boundary determination, and score calculation. More details of sub-network score calculations are described in Supplementary Note S5.

2.7 Experimental settings

2.7.1 Data collection and preprocessing

We assayed both tissue biopsy and liquid biopsy data in this study. First, for benchmark testing on cancer diagnosis and prognosis, we collected multiple datasets of different cancer types from TCGA (tissue data) to evaluate the classification performance, including 10 datasets for early- and late-stage classification, 10 datasets for low- and high-risk survival classification, and 5 datasets for drug responses prediction (Supplementary Fig. S2). In addition, we also collected and processed two types of body fluid datasets: the plasma dataset [373 samples assayed by total cell-free RNA-seq (ref. Chen , ref. Tao )] and the platelet dataset [918 samples assayed by blood platelet RNA-seq (ref. Best , ref. Best )]. Through our biological information pipeline, 3 and 7 biological modalities were derived from the TCGA (tissue biopsy) datasets and the liquid biopsy datasets, respectively. More details are described in Supplementary Notes S2.

2.7.2 Model training and test

We implemented Pathformer’s network architecture using the “PyTorch” package in Python v3.6.9, and our codes can be found in the GitHub repository (https://github.com/lulab/Pathformer). For model training and test, we used 5-fold cross-validation, and repeated it twice by shuffling. Before evaluating the performance on test sets, we optimized hyperparameters (e.g. learning rate, dropout probability of classification and constant coefficient for row-attention) and epoch numbers inside the training set only. More details of model training and test are described in Supplementary Note S6.

2.7.3 Evaluation criteria

When evaluating the classification performance, we used at least three evaluation indicators, area under the receiver operating characteristic curve (AUC), weighted-averaged F1 score (F1score_weighted), and macro-averaged F1 score (F1score_macro). Notably, we prioritized F1score_macro as the main evaluation criterion in this paper. This choice stems from the imbalance of sub-classes in our data, where F1score_macro stands out as a fairer and more robust indicator compared to other metrics such as AUC.

3 Results

3.1 Benchmark of Pathformer and 18 multi-omics data integration methods using TCGA data

We conducted a meticulous benchmark of Pathformer and 18 other multi-omics integration methods for various classification tasks in cancer diagnosis, treatment, and prognosis (Fig. 2). These methods can be categorized into three types. Type I includes early and late integration methods based on conventional classifiers, such as support vector machine (SVM), logistic regression (LR), random forest (RF), and extreme gradient boosting (XGBoost). Type II includes partial least squares-discriminant analysis (PLSDA) and sparse partial least squares-discriminant analysis (sPLSDA) of mixOmics (ref. Rohart ). Type III consists of deep learning-based integration methods, i.e. eiNN (ref. Preuer ), liNN (ref. Kuru ), eiCNN (ref. Fu ), liCNN (ref. Islam ), MOGONet (ref. Wang ), MOGAT (ref. Xing ), P-NET (ref. Elmarakeby ) and PathCNN (ref. Oh ). Among these, eiNN and eiCNN are early integration methods based on NN and CNN; liNN and liCNN are late integration methods based on fully connected neural network (NN) and convolutional neural network (CNN); MOGONet and MOGAT are multi-modal integration methods based on graph neural network; P-NET and PathCNN are representative multi-modal integration methods that combines pathway information. More details of comparison methods are in Supplementary Note S7.

To evaluate the performance, we tested the methods on multiple TCGA datasets for three tasks: cancer survival prediction, cancer staging, and drug response prediction. DNA methylation, DNA copy number variation (CNV), and RNA expression were used as input. Optimal hyperparameter combination for each dataset are listed in Supplementary Table S2. Considering the imbalanced numbers of sub-classes in the TCGA data, we utilized the macro-averaged F1 score as the primary evaluation metric for hyperparameter optimization and performance evaluation (Fig. 2 and Tables 1). Other evaluation indicators (e.g. AUC) are listed in Supplementary Table S3.

In general, Pathformer significantly performed better than 18 other integration methods in terms of F1score_macro score (Fig. 2a and b) and cross-validation variances (Fig. 2c). In cancer low- and high-risk survival classification tasks, comparing to the other eight deep learning methods (Type III), Pathformer’s F1score_marco showed average improvements between 6.3% and 14.6%. When comparing to eiXGBoost, which performed best in the conventional machine learning methods (Types I and II), Pathformer’s F1score_marco showed an average improvement of 8.3% (Fig. 2b). In early- and late-stage classification tasks, comparing to the deep learning methods (Type III), Pathformer’s F1score_marco showed average improvements between 5.1% and 12%. Compared to eiXGBoost, Pathformer’s F1score_marco showed an average improvement of 6.2% (Fig. 2b). In drug response prediction tasks, comparing to the deep learning methods (Type III), Pathformer’s F1score_marco showed average improvements between 8% and 13.6%. When comparing to eiXGBoost, Pathformer’s F1score_marco showed an average improvement of 8.6% (Fig. 2b). Moreover, Pathformer demonstrated reduced variance (Fig. 2c) and a stronger correlation between predictive confidence scores and fraction of positives (Supplementary Fig. S7) in cross-validation, indicating greater stability and reliability.

The detailed performance comparisons of Pathformer and other integration methods for different cancer types are shown in Tables 1 and Supplementary Table S3. In survival classifications, Pathformer achieved the highest F1score_macro and F1score_weighted in all the 10 datasets, and the highest AUC in 7 of 10 datasets. In stage classifications, Pathformer achieved the highest F1score_macro in 9 of 10 datasets, the highest F1score_weighted in 8 of 10 datasets, and the highest AUC in 6 of 10 datasets. In drug response prediction, Pathformer achieved the highest F1score_macro, F1score_weighted and AUC in all datasets.

3.2 Ablation analysis of Pathformer

We used ablation analysis to evaluate the contributions of different input modalities and calculation modules in the Pathformer model, based on nine datasets for cancer survival prediction (Fig. 3), nine datasets for cancer stage classification (Supplementary Fig. S8), and five datasets for drug response prediction (Supplementary Fig. S9). The pan-cancer dataset in cancer survival and stage classification was not used here. Firstly, to evaluate the contribution of integrating different modalities of data to classification, we compared seven models, including Pathformer with three modalities as input (RNA expression + DNA methylation + DNA CNV), Pathformer with two modalities as input (RNA expression + DNA methylation, RNA expression + DNA CNV, and DNA methylation + DNA CNV), and Pathformer with a single modality as input (RNA expression-only, DNA methylation-only, and DNA CNV-only). By comparing the performances of these models on cancer survival risk classification, we discovered that the model with all three modalities as input achieved the best performance, followed by the model with RNA expression and DNA CNV, and the model with DNA methylation-only (Fig. 3a). Furthermore, we observed that the performances of models with single modality as input can vary greatly between datasets. For example, DNA methylation-only model performed better than RNA expression-only and DNA CNV-only model in the LUSC, LIHC, and LGG datasets, but the opposite results were observed in the LUAD and BLCA datasets. Ablation analysis of different modalities on cancer stage classification (Supplementary Fig. S8a) and drug response prediction (Supplementary Fig. S9a) showed similar results. These findings underscore the distinct behaviors of different modalities in different cancer types, highlighting the necessity of multi-modal data integration in various cancer stage and survival risk classification tasks.

Next, to evaluate the essentialities of different calculation modules in Pathformer, we compared four additional variations of Pathformer, namely CC-attention, Transformer, PSNN, and NN, in which one to multiple modules of Pathformer are successively removed. The “CC-attention” model is Pathformer without pathway crosstalk network bias. The “Transformer” model is Pathformer without pathway crosstalk network bias and row-attention, using only normal attention mechanism and pathway embeddings. The “PSNN” model directly uses classification module with pathway embedding as input. The “NN” model directly uses classification module with gene embedding as input. As shown in Fig. 3b and Supplementary Figs S8 and S9, the complete Pathformer achieved the best classification performance, while the performance of CC-Attention, Transformer, PSNN, and NN decreased successively. This indicates that pathway crosstalk network, attention mechanism, and pathway embedding are all integral components of Pathformer. In particular, CC-attention exhibited significantly poorer classification performance compared to Pathformer, providing strong evidence for the necessity of incorporating pathway crosstalk in Pathformer.

3.3 Biological interpretability of Pathformer in breast cancer prognosis prediction using tissue data

To further understand the decision-making process of Pathformer and validate the reliability of its biological interpretability, we showed a case study on breast cancer survival risk classification. We demonstrated that Pathformer can use attention weights and SHAP values to identify modalities, pathways, and genes statistically associated with breast cancer prognosis, which aligns with known biological knowledge (Fig. 4).

First, at the omics and modality level, we visualized the contributions of different modalities for breast cancer survival risk classification by the attention weights (Fig. 4a). The contribution of transcriptomic data was the greatest in breast cancer prognostic prediction, which is consistent with the results of ablation analysis (Fig. 3a) and findings from other literatures (ref. Huang , ref. Tong ). Additionally, from Fig. 4a and Supplementary Fig. S10a, we observed that the contribution of various features in the same modality varied between BRCA prognosis and staging, such as DNA methylation. These findings further validate the necessity of biological multi-modal embedding and integration.

Next, at the pathway and gene levels, we identified key pathways with top 15 SHAP values and key genes with top 5 SHAP values for each pathway in breast cancer survival risk classification (Fig. 4b). Then, we presented a hub module of the updated pathway crosstalk network (Fig. 4c). These key pathways and genes identified by SHAP and the hub module are biologically meaningful and consistent with previous biological experiments. For instance, complex I biogenesis pathway, which was identified as the most critical pathway during the classification and a key node in the hub module of the updated pathway crosstalk network, was reported to play an important role in cancer cell proliferation and metastasis (ref. Urra ). Five mitochondrial genes (MT-ND4, MT-ND1, MT-ND3, MT-ND6, and MT-ND2), which were identified as key genes of the complex I biogenesis pathway by Pathformer, were also reported to be associated with breast cancer prognosis (ref. Kopinski ). Another example is FTL, which was predicted by our SHAP value to be up-regulated in the high-risk cancer group, was also reported to promote breast cancer cell proliferation validated by knockout experiments (ref. Tang ).

Subsequently, to facilitate a more intuitive understanding of the impact of active pathways identified by Pathformer on breast cancer survival risk classification, we depicted survival curves comparing patients with high and low scores in active pathways (Fig. 4d and Supplementary Fig. S11). The pathway score for each sample was obtained by averaging across different dimensions of pathway embedding updated by Pathformer. Log-rank tests indicated that most active pathways identified by Pathformer, like complex I biogenesis pathway, significantly influenced patient survival. Finally, to gain further insights into how Pathformer uses key features for accurate decision-making, we visualized pathway embedding changes and discussed the commonalities among correctly classified samples (Supplementary Note S9). We visually examined the feature extraction capability of Pathformer’s Transformer module with CC-attention by comparing PCA of pathway embedding matrices before and after Pathformer update, and pathway score heatmap (Supplementary Fig. S12a and b). We also found that Pathformer’s accuracy is not influenced by clinical indicators such as cancer subtype and age (Supplementary Fig. S12c).

3.4 Performance of Pathformer for the noninvasive diagnosis of cancer using liquid biopsy data

In clinical practice, cancer diagnosis involves not only using tissue data for cancer staging but also using liquid biopsy data (i.e. plasma) for noninvasive early detection and screening. The latter has even greater clinical significance because early detection substantially increases five-year survival rate of cancer patients. For instance, 5-year survival rates of colon cancer were reported as 93.2% for stage I, and only 8.1% for stage IV (ref. O’Connell ). Therefore, we applied Pathformer to liquid biopsy data, aiming to classify cancer patients from healthy controls. We curated two types of cell-free RNA sequencing (cfRNA-seq) data, including plasma datasets (comprising 98 healthy donors and 275 cancer samples) and platelet datasets (comprising 286 healthy donors and 632 cancer samples). We then calculated seven RNA-level modalities as Pathformer’s multi-modal input, including RNA expression, RNA splicing, RNA editing, RNA alternative promoter (RNA alt. promoter), RNA allele-specific expression (RNA ASE), RNA single nucleotide variations (RNA SNV), and chimeric RNA. Liquid biopsy data collection and preprocessing procedures are in Supplementary Note S2, while model parameters and settings are in Supplementary Note S6. Because these seven modalities of RNA may have information redundancy, we selected the best modality combination based on 2 times 5-fold cross validations (Supplementary Note S10). The results showed that the plasma data with seven modalities and the platelet data with three modalities obtained the best performances (AUCs > 0.9). Additionally, we found that Pathformer’s performance was superior to the other integration methods using the liquid biopsy data (Tables 2 and Supplementary Table S4). Because cancer screening usually requires high specificity, we particularly report sensitivities on 99% specificity in Table 2. Pathformer achieves an average sensitivity of 48.8% in the plasma dataset and an average sensitivity of 48.1% in the platelet dataset. It is worth noting that the sensitivity is still above 45% on 99% specificity in the plasma data even for the early-stage cancer patients, showing Pathformer’s potential for early cancer diagnosis.

3.5 Biological interpretability of Pathformer in the data of cancer patient’s blood

Based on the above analysis, we attempted to gain new insight into the deregulated alterations in plasma through Pathformer’s biological interpretability module (Fig. 5a and Supplementary Fig. S18a). First, we found that the pathways and genes ranked highly in SHAP values were associated with dysregulated alterations reported by previous experimental studies. For example, binding and uptake of ligands (e.g. oxidized low-density lipoprotein, oxLDL) by scavenger receptors pathway, with top SHAP value ranking, was reported to play a crucial role in cancer prognosis and carcinogenesis by promoting the degradation of harmful substances and accelerating the immune response (ref. Ryu ). Another two examples are DAP12 signaling pathway and DAP12 interactions pathway, which were highly ranked by SHAP value in both plasma and platelet data, were reported to regulate natural killer cell immune responses against certain tumor cells through platelet modulation cells (ref. Campbell and Colonna 1999, ref. Placke ).

Furthermore, Pathformer can explore potential novel interactions between various biological processes in cancer patients’ plasma by updating pathway crosstalk network (Fig. 5b). For example, the link between binding and uptake of ligands by scavenger receptors pathway and iron uptake and transport pathway was a novel addition to the known links (Pathformer’s input pathway crosstalk network curated from published databases, see Methods). This finding aligns with a previous report of SCARA5 (scavenger receptor class A member) as a ferritin receptor (ref. Yu ). The crosstalk between two pathways was amplified by Pathformer in plasma dataset, probably because they were important for classification. In summary, Pathformer’s updated pathway crosstalk network can effectively visualize the information flow between pathways related to cancer classification tasks, providing new insight into the crosstalk of biological pathways in cancer patients’ plasma.

4 Conclusion and discussion

Pathformer successfully applied a Transformer model to integrate multi-modal data for cancer diagnosis and prognosis. Particularly, it introduced a novel biological embedding method based on the compacted multi-modal vectors (Fig. 1b). Moreover, it utilized the criss-cross attention mechanism of Transformer to capture crosstalk between biological pathways and regulation between modalities (i.e. different omics).

4.1 Clinical applications of Pathformer

Pathformer can be applied to various classification tasks in disease diagnosis, treatment, and prognosis, such as early detection, cancer staging, survival prediction, and drug response prediction. Its predictive accuracy, stability, reliability, and biological interpretability were demonstrated through substantial benchmark and case studies focusing on cancer prognosis and noninvasive diagnosis. Ablation analysis demonstrated the pivotal role of various modal integrations and core modules (CC-attention, PNSS, and pathway embedding) in Pathformer for accurate classification. Our discussion on explained variances across integration models and multi-omics data further corroborated conclusions from benchmark tests and ablation analyses (Supplementary Note S11). Moreover, this framework is adaptable for the diagnosis and prognosis of other complex diseases, like autoimmune disease, neurodegenerative diseases, etc.

4.2 Potential targets revealed in cancer patients’ blood

Particularly, we identified some potential noninvasive cancer diagnostic biomarkers by Pathformer, such as the scavenger receptor related pathways and DAP12 related pathways, which are associated with extracellular vesicle transport (ref. Kzhyshkowska ) and immune response (ref. Campbell and Colonna 1999), respectively. We even found a new cancer-related pathway crosstalk in blood, which is between binding and uptake of ligands by scavenger receptors pathway and iron uptake and transport pathway. These results provide candidate targets for the mechanism study of cancer microenvironment and immune system, and even new targets for cancer treatment.

4.3 Limitations of Pathformer and future directions

For gene selection, Pathformer used genes involved in four common pathway databases, all of which consist of protein-coding genes. However, a substantial body of literature has reported that noncoding RNAs are also crucial in cancer prognosis and diagnosis (ref. Qi ). Therefore, incorporating noncoding RNAs and their related functional pathways into Pathformer would be promising for future work. For clinical applications in liquid biopsy, we used the multi-modal features derived from cfRNA-seq only in the application of liquid biopsy, because the published cell-free multi-omics datasets (ref. Tao ) are usually too small to be train-and-tested. For computational efficiency and memory costs, there is still room for improvement for Pathformer. Pathway embedding of Pathformer has prevented memory overflow of Transformer module caused by long inputs, but training still requires significant time and space (Supplementary Note S12). Therefore, when adding more pathways or gene sets (e.g. transcription factors), Pathformer still faces the issue of memory overflow. In the future work, we may introduce linear attention to further improve computational speed. Furthermore, potential signatures, regulations and biomarkers identified by Pathformer are also needed to be studied and validated by further biological experiments and clinical tests.

References

1. Swarm intelligence-enhanced detection of non-small-cell lung cancer using tumor-educated platelets.. Cancer Cell, 2017. [PubMed]
2. RNA-Seq of tumor-educated platelets enables blood-based pan-cancer, multiclass, and molecular pathway cancer diagnostics.. Cancer Cell, 2015. [PubMed]
3. DAP12: a key accessory protein for relaying signals by natural killer cell receptors.. Int J Biochem Cell Biol, 1999. [PubMed]
4. The cancer genome atlas pan-cancer analysis project.. Nat. Genet, 2013. [PubMed]
5. Cancer type classification using plasma cell-free RNAs derived from human and microbes.. Elife, 2022. [PubMed]
6. Predicting drug response of tumors from integrated genomic profiles by deep neural networks.. BMC Med Genomics, 2019. [PubMed]
7. Reactome: a database of reactions, pathways and biological processes.. Nucleic Acids Res, 2010. [PubMed]
8. Cui H , WangC, MaanH et al scGPT: toward building a foundation model for single-cell multi-omics using generative AI. Nature Methods 2024:1–11.
9. Biologically informed deep neural network for prostate cancer discovery.. Nature, 2021. [PubMed]
10. A gene prioritization method based on a swine multi-omics knowledgebase and a deep learning model.. Commun Biol, 2020. [PubMed]
11. PASNet: pathway-associated sparse deep neural network for prognosis prediction from high-throughput data.. BMC Bioinformatics, 2018. [PubMed]
12. Multi-omics approaches to disease.. Genome Biol, 2017. [PubMed]
13. SALMON: survival analysis learning with multi-omics neural networks on breast cancer.. Front Genet, 2019. [PubMed]
14. An integrative deep learning framework for classifying molecular subtypes of breast cancer.. Comput Struct Biotechnol J, 2020. [PubMed]
15. Highly accurate protein structure prediction with AlphaFold.. Nature, 2021. [PubMed]
16. KEGG: Kyoto encyclopedia of genes and genomes.. Nucleic Acids Res, 2000. [PubMed]
17. A hidden oncogenic positive feedback loop caused by crosstalk between Wnt and ERK pathways.. Oncogene, 2007. [PubMed]
18. Mitochondrial DNA variation and cancer.. Nat Rev Cancer, 2021. [PubMed]
19. MatchMaker: a deep learning framework for drug synergy prediction.. IEEE/ACM Trans Comput Biol Bioinform, 2022. [PubMed]
20. Stabilin‐1, a homeostatic scavenger receptor with multiple functions.. J Cell Mol Med, 2006. [PubMed]
21. A global pathway crosstalk network.. Bioinformatics, 2008. [PubMed]
22. Duality of interactions between TGF-β and TNF-α during tumor formation.. Front Immunol, 2021. [PubMed]
23. A unified approach to interpreting model predictions.. Adv Neural Inf Process Syst, 2017
24. A comprehensive evaluation of full-spectrum cell-free RNAs highlights cell-free RNA fragments for early-stage hepatocellular carcinoma detection.. EBioMedicine, 2023. [PubMed]
25. BioCarta.. Biotech Softw Internet Rep Comput Softw J Sci, 2001
26. Colon cancer survival rates with the new American Joint Committee on Cancer sixth edition staging.. J Natl Cancer Inst, 2004. [PubMed]
27. A novel method for crosstalk analysis of biological networks: improving accuracy of pathway annotation.. Nucleic Acids Res, 2017. [PubMed]
28. PathCNN: interpretable convolutional neural networks for survival prediction and pathway analysis applied to glioblastoma.. Bioinformatics, 2021. [PubMed]
29. Osseni MA , TossouP, LavioletteF, et al MOT: a multi-omics transformer for multiclass classification tumour types predictions. bioRxiv, 10.1101/2022.11.14.516459, 2022, preprint: not peer reviewed.
30. Modulation of natural killer cell anti-tumor reactivity by platelets.. J Innate Immun, 2011. [PubMed]
31. Opportunities and challenges provided by crosstalk between signalling pathways in cancer.. Oncogene, 2016. [PubMed]
32. DeepSynergy: predicting anti-cancer drug synergy with deep learning.. Bioinformatics, 2018. [PubMed]
33. Circulating long non-coding RNAs in cancer: current status and future perspectives.. Mol Cancer, 2016. [PubMed]
34. mixOmics: an R package for ‘omics feature selection and multiple data integration.. PLoS Comput Biol, 2017. [PubMed]
35. Scavenger receptor class a to E involved in various cancers.. Chonnam Med J, 2020. [PubMed]
36. PID: the pathway interaction database.. Nucleic Acids Res, 2009. [PubMed]
37. MOLI: multi-omics late integration with deep neural networks for drug response prediction.. Bioinformatics, 2019. [PubMed]
38. Overexpression of ferritin light chain as a poor prognostic factor for breast cancer.. Mol Biol Rep, 2023. [PubMed]
39. Cell-free multi-omics analysis reveals potential biomarkers in gastrointestinal cancer patients’ blood.. Cell Rep Med, 2023. [PubMed]
40. Undisclosed, unmet and neglected challenges in multi-omics studies.. Nat Comput Sci, 2021. [PubMed]
41. Transfer learning enables predictions in network biology.. Nature, 2023. [PubMed]
42. Integrating multi-omics data by learning modality invariant representations for improved prediction of overall survival of cancer.. Methods, 2021. [PubMed]
43. The mitochondrial complex (I) ty of cancer.. Front Oncol, 2017. [PubMed]
44. MOGONET integrates multi-omics data using graph convolutional networks allowing patient classification and biomarker identification.. Nat Commun, 2021. [PubMed]
45. Xing X , YangF, LiH, et al An interpretable multi-level enhanced graph attention network for disease diagnosis with gene expression data. In: 2021 IEEE International Conference on Bioinformatics and Biomedicine (BIBM). Houston, TX: IEEE, 2021, 556–61.
46. Interactions of ferritin with scavenger receptor class a members.. J Biol Chem, 2020. [PubMed]