Hybrid GNN-Transformer model for multi-omic cancer classification with interpretable pathway-driven feature selection

Data collection

The proposed dataset consists of gene data, methylation sites, and pathway activities. The gene data are collected from The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/). This includes expression values (RNA-seq), mutation status (binary: 1 if mutated, 0 otherwise), and copy number variations (CNV) status (e.g., −1 = deletion, 0 = neutral, +1 = amplification). The methylation sites are collected through UCSC Xena Browser but originating from TCGA, including beta-value (methylation level: 0 = unmethylated, 1 = fully methylated), genomic coordinates (chromosome, start, and end), and association with CpG islands/shores, with all measurements mapped to individual TCGA patients. The pathways data is curate cancer-related pathways collected from MSigDB database, representing biological relationships rather than patient-specific measurements. For cancer classification, the dataset is composed of 6,356 samples of 23 cancer types, with 80% (5,085 samples) for training, 10% (635 samples) for validation, and 10% (636 samples) for testing. For breast cancer subtypes classification, the training dataset is composed of 974 samples of five breast cancer subtypes, and the test dataset is composed of 122 samples collected from the Clinical Proteomic Tumor Analysis Consortium (CPTAC).

Data preprocessing

The collected data is preprocessed in order to ensure effective learning of the proposed model. The gene expressions are normalized using log normalization.

(1) $$G_{log}={\log }_2(G+1).$$The resulted expressions are further normalized using z-score normalization,

(2) $$G_{zscore}=\frac{(G_{log}-{\mu }_G)}{{\sigma }_G},$$where ${\mu }_G$ and ${\sigma }_G$ are the mean and the standard deviation over all samples, respectively. The beta values of the DNA methylation data are first converted to M-values.

(3) $$M_{values}={\log }_2\left(\frac{\beta }{1-\beta +\epsilon }\right).$$

High-dimensional DNA methylation profiles, with $\sim 450$ k CpG sites, presents significant computational challenges for cancer classification. Traditional dimensionality reduction techniques like Principal Component Analysis (PCA) often fail to preserve critical spatial patterns in methylation data (e.g., CpG island clusters), which are biologically significant.

In this article, we propose a hybrid approach combining HWT with PSO to achieve efficient feature reduction while retaining diagnostically relevant information. HWT is a multi-scale feature extraction method that decomposes signals into hierarchical approximations (low-frequency trends) and details (high-frequency variations). HWT can decompose data into multi-resolution coefficients, capturing both local and global patterns. In general, HWT is applied iteratively at multiple scales, that is, at each scale, the signal is downsampled, and the resulting coefficients become the input for the next scale. Thus, each scale captures different levels of detail in the signal. In this work, HWT is utilized to decompose methylation M-values (treated as a 1D signal) into approximation $A_j$ (low-frequency trends or broad methylation domains) and detail $D_j$ (high-frequency noise or fine-scale variations) coefficients, where $j$ represents the level of the wavelet decomposition. A three-level decomposition of the signal ( $j=1,2,$ and $3$) is carried out with high-resolution details (individual CpG sites), mid-resolution details (gene clusters), and low-resolution approximation (chromosomal regions). For a methylation profile $M=[m_1,m_2,\cdots ,m_N]$, HWT is computed by

(4) $$A_1=[\frac{m_1+m_2}{\sqrt{2}},\frac{m_3+m_4}{\sqrt{2}},\cdots ,\frac{m_{N-1}+m_N}{\sqrt{2}}],$$

(5) $$D_1=[\frac{m_1-m_2}{\sqrt{2}},\frac{m_3-m_4}{\sqrt{2}},\cdots ,\frac{m_{N-1}-m_N}{\sqrt{2}}],$$and

(6) $$A_j=[\frac{A_{j-1,1}+A_{j-1,2}}{\sqrt{2}},\frac{A_{j-1,3}+A_{j-1,4}}{\sqrt{2}},\cdots ,\frac{A_{j-1,N/2^{j-1}-1}+A_{j-1,N/2^{j-1}}}{\sqrt{2}}],$$

(7) $$D_j=[\frac{A_{j-1,1}-A_{j-1,2}}{\sqrt{2}},\frac{A_{j-1,3}-A_{j-1,4}}{\sqrt{2}},\cdots ,\frac{A_{j-1,N/2^{j-1}-1}-A_{j-1,N/2^{j-1}}}{\sqrt{2}}],$$where $j=2$ and $3$. Thus, high-frequency detail coefficients, $D_j$ (for levels $j=1,2,$ and 3) are discarded to eliminate noise, while approximation coefficients $A_j$ are progressively compressed, maintaining spatial genomic integrity. Through three-level compression, $A_1$ contains N/2 coefficients, $A_2$ contains $N/4$ coefficients, and $A_3$ retains only $N/8$ coefficients (i.e., 12.5% of N). As a result, the utilized HWT decomposition reduces the features by 87.5% (i.e., from 450 to 56 K coefficients).

While HWT reduces dimensionality, not all wavelet coefficients are discriminative for cancer subtypes. We combined HWT with PSO to adaptively identify and retain maximally discriminative coefficients. PSO is a metaheuristic optimization technique inspired by swarm intelligence, where candidate solutions (“particles”) navigate high-dimensional search spaces. Each particle iteratively refines its position based on personal experience ( $p^{best}$) and swarm knowledge ( $g^{best}$), dynamically balancing exploration of new solutions and exploitation of known successes. This enables efficient convergence toward coefficient subsets that maximize inter-class separation between cancer subtypes. We employ PSO to identify the most informative subset of HWT coefficients. Specifically, each particle encodes a binary vector ${\mathbf{P}}_i\in \{0,1{\}}^{N_w}$, where $N_w$ is the number of wavelet coefficients, where a value of 1 indicates the coefficient is selected. Each particle’s feature subset is evaluated using a support vector machine (SVM) classifier with a radial basis function (RBF) kernel trained solely on the training subset. The classification accuracy, computed as the average accuracy from three-fold cross-validation on the training data, is then used as a fitness function, whose input is element-wise multiplication of the wavelet-transformed data and ${\mathbf{P}}_i$, to determine each particle’s personal best solution ( $p_{i,d}^{best}$) and the swarm’s global best solution ( $g_d^{best}$). Based on these individual and swarm best solutions, the particles adjust their positions during optimization:

(8) $$v_{i,d}^{t+1}=w{v}_{i,d}^t+c_1{r}_1(p_{i,d}^{best}-p_{i,d}^t)+c_2{r}_2(g_d^{best}-p_{i,d}^t),$$

(9) $$p_{i,d}^{t+1}=\{\begin{array}{cc}1\hfill & if\sigma (v_{i,d}^{t+1})>0.5\hfill \\ 0\hfill & otherwise\hfill \end{array}$$where $\sigma$ is the Sigmoid function, $w$ is inertia, $c_1$ and $c_2$ are learning rates, $r_1$, and $r_2$ are random numbers. Note that $i$ denotes the particle index, $d$ denotes the index of the binary vector ${\mathbf{P}}_i$, and $t$ denotes the iteration number. Thus, PSO is used to optimize the features subset that maximizes cancer subtype separability. To prevent data leakage, the dataset was divided into training (80%), validation (10%), and test (10%) subsets prior to any feature selection. The PSO feature selection process operated exclusively on the training subset, with three-fold cross-validation performed within this subset. The optimal features obtained from PSO were then applied to all subsets, i.e., training, validation, and test subsets.

Heterogeneous graph representation

To capture the complex interactions underlying cancer biology, we aim to build a graph that represents genes, methylation sites, and pathways as nodes, with edges encoding their biological relationships (e.g., regulation, membership, interaction). That is, each node has its own type and attributes as given below.

• -

  Genes: This node represent genes, which are functional units of DNA that encode proteins or regulatory RNAs. The node attributes are expression levels (RNA-seq), mutation status, and CNV status. For example, the gene “TP53” has the attributes, expression: 12.5 TPM (high in normal tissue, low in tumor), mutation: 1 (mutated in 50% of cancers), and CNV: 0.

• -

  Methylation sites: This node represents CpG sites, which are regions of DNA where a cytosine nucleotide is followed by a guanine nucleotide. Methylation at these sites regulates gene expression. The node attributes are beta-values, genomic coordinates (chromosome, start, end), and association with CpG islands/shores. For example: the CpG site “cg03607383” has the features, beta-value: 0.9 (hypermethylated), and location: promoter of BRCA1 (chr17:43,044,294).

• -

  Pathways: This node represents biological pathways, which are collections of genes that work together to perform specific functions (e.g., cell cycle regulation). The node attributes are pathway size (number of genes), functional category (e.g., “Metabolism,” “Immune Response”), and activity score (computed via AUCell). For example: the pathway “KEGG”, hsa04110 (Cell Cycle), genes: “CDK1”, “CCND1”, “TP53”, and activity score: 2.3 (upregulated in tumors).

The edges of the proposed heterogeneous represent biological relationships between nodes, with each edge type encoding a specific biological interaction. The description of these edges is given below.

• -

  Gene-gene: It represents physical or functional interactions between genes (e.g., protein-protein interactions, PPI, and co-expression). For PPI, only high-confidence interactions (STRING score > $0.7$) are used. For example, “TP53” $\to$ “MDM2” (“MDM2” ubiquitinates “TP53” for degradation). For co-expression, only pearson correlation > $0.6$ (or $<-0.6$ for negative correlation) is used. For example, the genes“EGFR” and “KRAS” are co-expressed in lung adenocarcinoma.

• -

  Gene-pathway: It represents the membership of genes in biological pathways. For example, the gene “TP53” is a key gene in the pathway “hsa04115 (p53 signaling)”.

• -

  Methylation-gene: it represents the regulatory effect of methylation on gene expression. For example, hypermethylation of the CpG site “cg03607383” silences “BRCA1” expression.

• -

  Methylation-pathway: it represents the indirect effect of methylation on pathways via target genes. That is, if a CpG site regulates a gene in a pathway, connect the site to that pathway. For example, hypermethylation of the CpG site “cg03607383” silences the gene “BRCA1”, disrupting the “DNA repair” pathway.

Thus, the heterogeneous graph is defined by

(10) $$G_i=(V,E),$$where

(11) $$V=gene\mathrm{\_}nodes\cup meth\mathrm{\_}nodes\cup pathway\mathrm{\_}nodes,$$and

(12) $$E=gene-gene,gene-pathway,meth-gene,meth-pathwayedges.$$

An example of heterogeneous subgraph is illustated in Fig. 2. Note that the blue, green, and red circles are corresponding to gene, methylation site, and pathway, respectively.

Graph neural network (GNN)

The GNN operates on the heterogeneous graph to capture local interactions between nodes, such as gene-gene interactions, methylation-gene regulation, and gene-pathway memberships. Rational graph convolutional network (RGCN) is an extension of the standard graph convolutional network (GCN) that handles heterogeneous graphs by using separate weight matrices for each edge type. This makes it ideal for modeling multi-omic data, where nodes and edges represent different biological entities and relationships. The core idea of RGCN is to aggregate information from neighbors while considering the type of relationship between nodes. For node $v$ and edge type $r$, the message-passing rule is:

$$h_v^{(l+1)}=\sigma \left(\sum _{r\in R}\sum _{u\in {\mathcal{N}}_r(v)}\frac{1}{|{\mathcal{N}}_r(v)|}{W}_r^{(l)}{h}_u^{(l)}+W_0^{(l)}{h}_v^{(l)}\right)$$where $h_v^{(l)}$ is embedding of node $v$ at layer $l$, ${\mathcal{N}}_r(v)$ are neighbors of $v$ via edge type $r$, $W_r^{(l)}$ are learnable weights for edge type $r$, and $\sigma$ is activation function.

RGCN in the proposed model consists of two layers to balance model complexity and computational efficiency. The raw node features are fed into the first layer to aggregate information from neighbors of each edge type and update node embeddings using the message-passing rule. The output of the first layer involves intermediate node embeddings capturing local interactions, which are fed into the second layer to perform another round of aggregation and updating. The output of the second layer is final node embeddings capturing higher-order interactions (e.g., pathways influenced by multiple genes).

In RGCN, each edge type has its own weight matrix, allowing the model to learn relationship-specific patterns. RGCN include self-loop connections to preserve node-specific information. The benefits of RGCN come from its seamless handling multiple node and edge types, its efficient processing of large graphs with sparse connections, and its ability to learn relationship-specific embeddings, making it easier to interpret results.

Transformer for context-aware feature selection

The transformer encoder processes the preprocessed data to capture global dependencies using pathway-guided attention. Gene expression and methylation sites are jointly input to the transformer to model their direct regulatory relationships (e.g., promoter methylation suppressing gene expression), enabling self-attention to capture cross-modal dependencies critical for cancer biology. The input is formed by concatenating the preprocessed genes and methylation sites metrices into a single matrix:

$$X=[X_{gene};X_{methylation}]\in {\mathbb{R}}^{N\times (D_g+D_m)}$$where N is number of samples, $D_g$ is number of selected gene features, and $D_m$ is number of selected methylation features. The preprocessed data X is projected into a shared feature space using learned embeddings:

(13) $$E=X{W}_E+b_E,$$where $W_E\in {\mathbb{R}}^{(D_g+D_m)\times D}$ is learnable embedding matrix, $b_E\in {\mathbb{R}}^D$ is learnable bias term, and D is dimension of the transformer embeddings. The architecture of the transformer encoder consists of the following components:

• Positional encoding: Since transformers are permutation-invariant, positional encodings are added to the input embeddings to preserve the order of features: (14) $$E_{pos}=E+P,$$where $P\in R^{(D_g+D_m)\times D}$ is a positional encoding matrix. Here, tokens are ordered by genomic coordinates, grouped by chromosome and sorted by ascending position, with genes and methylation sites interleaved to preserve spatial regulatory relationships. Physically proximate elements become consecutive tokens, receiving similar sinusoidal positional encodings that model cis-regulatory interactions.

• Self-attention mechanism: This mechanism computes attention scores between all pairs of features, allowing the model to focus on the most relevant features for each prediction. The input embeddings $E_{pos}$ are projected into query (Q), key (K), and value (V) matrices: $$Q=E_{pos}{W}_Q,K=E_{pos}{W}_K,V=E_{pos}{W}_V,$$where $W_Q,W_K,W_V\in {\mathbb{R}}^{D\times D_k}$ are learnable weight matrices. The attention scores are computed by $$\mathrm{A}\mathrm{t}\mathrm{t}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}(Q,K,V)=\mathrm{s}\mathrm{o}\mathrm{f}\mathrm{t}\mathrm{m}\mathrm{a}\mathrm{x}\left(\frac{Q{K}^T}{\sqrt{d_k}}\right)V,$$where $d_k$ is dimension of the key vectors and used for scaling, and $softmax$ normalizes the scores to sum to 1.

• Pathway-guided masking: To ensure biological relevance, a binary mask M is applied to the attention matrix as $$Attention\mathrm{\_}weights=Attention\odot M,$$where $M_{i,j}=1$ if features $i$ and $j$ are biologically related (e.g., genes in the same pathway), and 0 otherwise. This mask is computed from RCGN pathway-informed embeddings: (15) $$M_{ij}=\{\begin{array}{cc}1\hfill & if\frac{{\displaystyle h_i\cdot h_j}}{\Vert h_i\Vert \Vert h_j\Vert }>\tau \vee i=j\hfill \\ 0\hfill & otherwise\hfill \end{array}$$where $h_i$ and $h_j$ are RCGN embeddings corresponding to positions in transformer input, and $\tau$ is a threshold dynamically computed per sample as the 75th percentile of all pairwise cosine similarities between RCGN embeddings in the transformer input, $\tau =\mathrm{q}\mathrm{u}\mathrm{a}\mathrm{n}\mathrm{t}\mathrm{i}\mathrm{l}{\mathrm{e}}_{0.75}(\cos (h_i,h_j)\mathrm{\forall }i\ne j)$.

• Multi-head attention The transformer uses multi-head attention to capture different types of relationships as $$MultiHead(Q,K,V)=CONCAT(hea{d}_1,hea{d}_2,hea{d}_H)W_o,$$where H is number of attention heads, $hea{d}_i=Attention(Q{W}_Q^i,K{W}_K^i,V{W}_V^i)$, and $W_o\in {\mathbb{R}}^{H{D}_k\times D}$ is a learnable projection matrix. Following this, layer normalization and residual connections are used: $$h_{attn}=LayerNorm(X+MultiHead(Q,K,V)).$$

• Feed-Forward Network (FFN) After multi-head attention, the output is passed through a position-wise feed-forward network as $$h_{ffn}=ReLU(h_{attn}{W}_1+b_1)W_2+b_2,$$where $W_1,W_2\in {\mathbb{R}}^{D\times D}$ and $b_1,b_2\in {\mathbb{R}}^D$ are learnable parameters. This is followed by layer normalization and residual connections to stabilize training as $$h_{transformer}=LayerNorm(h_{ffn}+h_{attn}).$$

In our model, the transformer encoder uses pathway-guided attention to ensure that attention is focused on biologically relevant relationships. This is achieved by restricting attention to genes and methylation sites within the same pathway and highlighting important features based on pathway membership. The RGCN and transformer modules are trained simultaneously in an end-to-end framework rather than as isolated components. This joint optimization enables dynamic interplay between local graph-based feature learning (RGCN) and global context modeling (transformer), where gradient signals from the classification loss co-adapt both modules. Specifically, the RGCN’s pathway-informed embeddings constrain the transformer’s attention mechanism via biological priors, while the transformer’s context features provide feedback to refine graph representations. This co-training strategy preserves biological consistency between local pathway structures and global feature interactions, which is critical for capturing cancer-specific dysregulation patterns. The local (RGCN) and global (transformer) features are fused for final predictions,

$$h_{fused}=CONCAT(h_{RGCN},h_{transformer})$$where $h_{RGCN}$ are node embeddings from the RGCN and $h_{transformer}$ are context-aware embeddings from the transformer.

Classification head

The classification head is the final component of the proposed model, which is responsible for predicting cancer subtypes based on the fused representation of local (RGCN) and global (transformer) features. It consists of a multi-layer perceptron (MLP) that maps the fused features to the predicted cancer subtypes.

$$\hat{y}=Softmax(W_{2}ReLU(W_1{h}_{fused}+b_1)+b_2).$$The cross-entropy loss is used for training, which is given by

$$L=\sum _i^N{y}_i\log ({\hat{y}}_i),$$where $y_i$ is the true label and ${\hat{y}}_i$ is the predicted probability for the $i$th sample.

Finally, it should be pointed out that the dual use of pathway data, in graph construction and pathway-guided attention, is not redundant, but a deliberate design to enforce biological plausibility at multiple levels. In graph construction, pathways define the topology of the heterogeneous graph. They establish hard constraints (e.g., gene-pathway memberships) that ensure the RGCN learns biologically grounded local interactions. Without this, the graph could capture spurious correlations. In pathway-guided attention, they act as a soft regularizer for global attention. While the graph uses pathways to define edges, the attention mechanism uses them to dynamically weight gene-gene interactions in the Transformer. Crucially, it masks non-pathway-related interactions (e.g., between genes in unrelated pathways) that might otherwise dominate due to statistical biases.

The pseudo-code for the proposed hybrid graph-transformer model for cancer classification is given in Algorithm 1.

Performance evaluation

The performance of the proposed method is evaluated on two classification tasks. The first task involves cancer classification using the data of 6,356 samples with 23 cancer types. The second task involves breast cancer subtypes classification, i.e., luminal A, luminal B, HER2-enriched, basal-, and normal-like (Parker et al., 2009). The training dataset consists of a total of 974 samples of breast cancer. A breast cancer cohort of 122 samples collected from CPTAC is used as a independent test dataset. For each class, the classification outcomes are categorized into true positive (TP) indicating samples that correctly diagnosed as belonging to this class, true negative (TN) indicating samples that correctly diagnosed as not belonging to this class, false positive (FP) indicating samples that wrongly diagnosed as belonging to this class, and false negative (FN) indicating samples that wrongly diagnosed as not belonging to this class. The evaluation is then carried out for each class independently using the performance metrics,

(16) $$\mathrm{A}\mathrm{c}\mathrm{c}\mathrm{u}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{y}(\mathrm{A}\mathrm{C}\mathrm{C})=\frac{TP+TN}{TP+TN+FP+FN},$$

(17) $$\mathrm{P}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{i}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}(\mathrm{p}\mathrm{o}\mathrm{s}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e}\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}\mathrm{i}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e}\mathrm{v}\mathrm{a}\mathrm{l}\mathrm{u}\mathrm{e},\mathrm{P}\mathrm{P}\mathrm{V})=\frac{TP}{TP+FP},$$

(18) $$\mathrm{S}\mathrm{p}\mathrm{e}\mathrm{c}\mathrm{i}\mathrm{f}\mathrm{i}\mathrm{c}\mathrm{i}\mathrm{t}\mathrm{y}(\mathrm{S}\mathrm{P}\mathrm{C})=\frac{TN}{TN+FP},$$

(19) $$\mathrm{S}\mathrm{e}\mathrm{n}\mathrm{s}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{i}\mathrm{t}\mathrm{y}(\mathrm{S}\mathrm{E}\mathrm{N})=\frac{TP}{TP+FN},$$

(20) $$\mathrm{F}1\text{–}\mathrm{s}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}=\frac{2\cdot \mathrm{P}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{i}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}\cdot \mathrm{S}\mathrm{e}\mathrm{n}\mathrm{s}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{i}\mathrm{t}\mathrm{y}}{\mathrm{P}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{i}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}+\mathrm{S}\mathrm{e}\mathrm{n}\mathrm{s}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{i}\mathrm{t}\mathrm{y}}.$$

Afterwards, the averages of these metrics are calculated as indicators for the overall prediction performance.

Classification of cancer types

First, the performance of the proposed model is evaluated on cancer classification task with respect to several baselines in the literature. These baselines are given below.

• -

  MOGTN (Kaczmarek et al., 2021): This method integrates miRNA and mRNA expression data into a biological interaction graph, where nodes represent molecular entities (miRNAs/mRNAs) and edges encode regulatory relationships derived from TargetScan databases.

• -

  MOGONET (Wang et al., 2021): This method utilizes GCN for omics-specific learning and unifies it with multi-omics integrative classification at the label space.

• -

  DeePathNet (Cai et al., 2024): A transformer-based explainable deep learning model integrating cancer-specific pathway information into multi-omic data analysis.

• -

  MOGNN (Li & Nabavi, 2024): An end-to-end multi-omics GNN framework utilizing multi-omics data in the form of heterogeneous multi-layer graphs, which combine both inter-omics and intra-omic connections from established biological knowledge. This model incorporates learned graph features and global genome features for accurate classification.

For fair comparison, all models were trained and evaluated on the dataset of 6356 TCGA samples across 23 cancer types, with strictly matched data partitioning: 80% (5,085 samples) for training, 10% (635 samples) for validation, and 10% (636 samples) for testing. However, modality inputs were adapted to each model’s capabilities as shown in Table 1, with our framework uniquely leveraging the complete multimodal feature set (gene expressions, methylation sites, and pathway activities). Stratified sampling preserved class distributions across splits.

The results are presented in Table 2, demonstrating that our model consistently outperforms the others in terms of all performance metrics. Further, the confusion matrix for our model is shown in Fig. 3 to illustrate the performance for each cancer type. As can be seen, our model achieves an accuracy higher than 97% for most cancer types.

Notably, breast invasive carcinoma (BRCA), acute myeloid leukemia (LAML), ovarian serous cystadenocarcinoma (OV), pheochromocytoma paraganglioma (PCPG), prostate adenocarcinoma (PRAD), skin cutaneous melanoma (SKCM), thyroid carcinoma (THCA), and uterine carcinsarcoma (UCS) are the most accurately classified cancer types. In contrast, rectum adenocarcinoma (READ) is the cancer type with the poorest accuracy, where 34% of the samples are misclassified as colon adenocarcinoma (COAD). This is owing to the fact that the colon and rectum are adjacent tissue types sharing highly similar attributes, where they often grouped together and are treated with similar chemotherapeutic regimens (The Cancer Genome Atlas Network, 2012). Also, the accuracy is stomach adenocarcinomas (STAD) is poor, where 13% of its samples misclassified as esophageal carcinoma (ESCA). This can be attributed to the similar histopathology and the anatomical proximity of STAD and ESCA (Akiyama et al., 1997).

We conducted a comprehensive ablation study to evaluate the contributions of each component in our hybrid graph-transformer framework, testing nine distinct model variants that systematically removed or altered key elements: pathway attention mechanisms, RGCN module for local interactions, transformer module for global context, gene expression data, methylation sites, pathway database, and HWT-PSO. Each variant was rigorously evaluated using standard classification metrics (Accuracy, F1-score, and AUC), a specialized BioConsistency metric (0–1 scale) quantifying alignment between model attention weights and established biological pathways, and training efficiency measured by epochs to convergence.

The results shown in Table 3 reveal critical dependencies and performance impacts across components. Removing the pathway attention mechanism caused a substantial 3.25% accuracy decline and severe 0.17 BioConsistency reduction, demonstrating its essential role in enforcing biologically plausible feature interactions. Eliminating the RGCN module resulted in the most dramatic performance degradation (8.12% accuracy drop), confirming its irreplaceable function in capturing local molecular interactions and pathway-constrained relationships. Disabling the transformer module reduced accuracy by 5.64%, highlighting its unique capacity to model global gene dependencies across tissue types. Data modality ablation showed that gene expression contributes marginally more predictive power than methylation (3.53% vs. 2.82% accuracy drop when removed), but both are necessary for comprehensive regulatory modeling. Strikingly, removing the external pathway database caused the second-largest accuracy reduction (7.04%) and most severe BioConsistency collapse (0.27 decrease), underscoring how biological knowledge grounds the entire architecture in reality. The 3.76% accuracy drop when using HWT without PSO confirms that wavelet decomposition alone preserves spatial patterns but fails to select biologically relevant features. PSO is essential for identifying discriminative wavelet coefficients that optimize cancer type separation.

Breast cancer subtypes classification

As a use case, we evaluate the performance of the proposed model against the baselines in classifying breast cancer subtypes, luminal A, luminal B, HER2-enriched, basal-, and normal-like (Parker et al., 2009). The training of the compared models was performed using a total of 974 samples of breast cancer. A breast cancer cohort of 122 samples collected from CPTAC is used as a independent test dataset. The results are shown in Table 4. It can be seen that our model achieves the best performance metrics in both training and test datasets. Further, the confusion matrix of the proposed model is shown in Fig. 4 to evaluate the performance for each subtype. Apparently, our model has its highest accuracy in the classification of the “Basal-like” subtype (98.5%). In contrast, the “Normal-like” subtype is poorly classified, where half of its samples are misclassfied as “Luminal A”. This is unsurprising due to the fact that both “Luminal A” and “Normal-like” share similar IHC markers (Cai et al., 2024); indeed, the “Normal-like” subtype is rarely used in medical sector (Raj-Kumar et al., 2019).

Enrichment analysis

To validate the biological and clinical significance of biomarkers prioritized by our hybrid GNN-Transformer model, we perform enrichment analysis on genes, methylation sites, and pathways associated with the basal-like breast cancer subtype. Pathway importance scores and SHAP values are used to identify drivers of basal-like pathogenesis and therapeutic vulnerabilities. Pathway importance scores highlight systems-level drivers, where high-score pathway nominating therapeutic targets and low-score pathways exclude irrelevant mechanisms. SHAP attributes the prediction to all features, assigning each feature an importance value.

The outcomes of the proposed model are shown in Fig. 5, where pathway importance scores are depicted in the top panel and SHAP values are depicted in the bottom panel. For pathway importance scores, our model prioritizes pathways with known roles in basal-like cancer biology. That is, dysfunctional DNA repair, particularly defective homologous recombination due to BRCA1 mutations, plays a key role in basal-like breast cancer. PI3K-AKT, which regulates cell survival, metabolism, and growth, is often activated in basal-like breast cancer due to mutations or amplifications in the genes encoding components of this pathway. The p53 signaling pathway is critically disrupted in basal-like breast cancer, with profound implications for tumor aggressiveness, genomic instability, and therapeutic resistance. For SHAP values, our model recapitulates known basal-like drivers (e.g., TP53 and CCNE1) and nominates novel candidates (e.g., PD-L1 promoter methylation). In other words, SHAP values highlight oncogenic drivers (CCNE1 and MYC) and epigenetic regulators (e.g., BRCA1 promoter methylation) as critical basal-like features. These findings align with aggressive phenotype and therapeutic resistance of basal-like cancer. As can be seen, only positive values are shown on the x-axis, representing the overall influence of each feature rather than its directional effect. This allows for more reliable identification of consistently impactful features, which is essential in the context of biomarker discovery.