paper.md

Abstract

Introduction:

Encoding prior biological knowledge directly into neural network architecture has emerged as a powerful strategy for prediction from high-dimensional molecular data. Rather than learning from a fully connected parameterization, biologically informed neural networks (BINNs), constrain inter-layer connectivity to follow curated gene–pathway relationships. Concretely, each the weight matrix is reparamatized where the binary mask encodes gene-to-pathway ancestry drawn from resources such as Reactome. Because every hidden unit corresponds to a named biological entity (gene, process, pathway, etc) , the resulting models are interpretable by construction.

A structural limitation of these models is that pathway-masked networks process each observation independently but never between them. Single-cell measurements however are sparse and noisy, dominated by dropout arising from low per-cell sequencing coverage and substantial biological variation. The standard response across the field is to assume that cells occupy a low-dimensional manifold, an assumption underlying dimensionality-reduction and clustering pipeline in routine use (reference them as well).

A Graph Neural Network (GNN), performs this step inside the network rather than a preprocessing step, attenuating per-cell noise while preserving shared signal.
Whether such averaging aids prediction depends on the relationship between the graph structure and the outcome label. Most GNN’s rely on the homophily assumption and heterophilous graphs, in which dissimilar nodes are linked, are widely held to limit the performance of architectures designed under that assumption [10]. Graph convolutions can perform well under heterophily provided different classes retain distinguishable neighborhood distributions, and degrade principally when the graph carries no class-discriminative information at all [9]. We therefore report the label homophily of each cell graph alongside predictive performance, rather than presupposing it.

We address the question of GNN’s performance with a minimal comparison constructing pathway-masked networks that are identical in biological priors and differing only in the presence of a normalized propagation operator over the cell-similarity graph (Methods). The identity as operator recovers standard BINN[1] with an adversarial neural network (ANN) encoder whereas retaining it yields a graph-convolutional (GCN) or graph-attention (GAT) variant. The three architectures are valuated on single-cell, bulk and spatial datasets under three trials of five-fold stratified cross-validation.

A common critique of GNN performance concerns transductive training, where edges between training and testing nodes allow information to leak from the test set during the graph's propagation step. To address this, we evaluate the models under both the full graph (transductive learning) and fold-dependent masking, which prevents propagation from out-of-set observations. For interpretability we apply tools such as integrated gradients and SHAP, and find no significant difference from traditional BINNs. The architectures, masks, and training loops are released as a general-purpose PyTorch library built on PyTorch Lightning's compute infrastructure, available at https://github.com/CompBio-Lab/BINN.

Results:

On the Image-based spatial transcriptomics dataset GSE250346[11] (n = 330,991 cells, 2 classes, positive-class prevalence 63.5%, edge homophily h = 0.9998), the BINN baseline (pathway-pruned ANN, no graph operator) reached AUROC 0.530 ± 0.050, versus 0.920 ± 0.040 (GCN) and 0.946 ± 0.037 (GAT with LayerNorm). The graph operator produced a large and consistent gain over the baseline (GCN: Δ = +0.390, paired Wilcoxon signed-rank test; GAT-LN: Δ = +0.416, p = 0.063), while the difference between the two graph variants was small (Δ = +0.026, p = 0.063 Wilcoxon / p = 0.002 paired t). Notably, GAT with BatchNorm collapsed to AUROC 0.675 ± 0.241 with 3 of 5 folds degenerating to single-class prediction, ranking below the simpler GCN; LayerNorm eliminated this collapse entirely (0/15 folds across 3 seeds), [TODO: provide explanation for collapse prevention.] The near-perfect homophily (h = 0.9998, only 299 of 1,985,946 edges crossing the class boundary) means the spatial graph is effectively carrying label information: because pathology annotations are assigned to spatially contiguous polygons, a cell’s KNN neighbors almost always share its class, so message passing functions as near-perfect label smoothing — which explains why both graph operators dominate the non-graph baseline but barely differ from each other.

On the bulk PAAD dataset Mulgo[] with 175 patients and 2 classes, the BINN baseline reached AUROC 0.697 ± 0.091, versus 0.800 ± 0.051 (GCN) and 0.806 ± 0.058 (GAT). The graph operator produced a consistent and statistically significant gain over the baseline, while the difference between the two graph variants was negligible (Δ = +0.006, 95% CI [−0.028, +0.039], paired t p = 0.73). The corresponding AUPR gains echo this pattern with a trial mean improvement from 0.673 to 0.778 and 0.775 for GCN and GAT respectively. The graph's moderate homophily h = 0.6408 means the cosine-kNN neighborhoods are meaningfully but far from perfectly aligned with the recurrence label; message passing therefore acts as a weak; informative label smoothing rather than the near-perfect propagation seen in the spatial-transcriptomics setting.

Evaluate on the CITE-seq dataset of 174,753 cells from 41 donors, the advantage of Graph BINN’s remain. TODO: will write

Add multi figure plot of results

References:

1. Biologically informed deep neural network for prostate cancer discovery (P-NET)

2. MULGONET: An interpretable neural network framework to integrate multi-omics data for cancer recurrence prediction and biomarker discovery

3. Interpreting biologically informed neural networks for enhanced proteomic biomarker discovery and pathway analysis (BINN)

4. A systematic review of biologically-informed deep learning models for cancer,

5. Applying interpretable machine learning in computational biology — pitfalls, recommendations and opportunities,

6. Graph Attention Networks Velickovic

7. Reactome: a database of reactions, pathways and biological processes

8. Gene Ontology: tool for the unification of biology

9. Is Homophily a necessity for graph neural networks?

10. Graph Neural Networks for Graphs with Heterophily: A Survey.

11. GSE250346 — Vannan, Lyu et al., Nature Genetics 57, 647–658 (2025).

12. A Critique on Transductive Evaluation for GNN Node Classification

13. GraphSCI: A graph neural network based imputation method for single-cell RNA-seq data. Genome Biology.

Methods:

Preprocessing

The spatial transcriptomics cohort comprises 330,991 cells from 28 patients of the GSE250346 Xenium pulmonary fibrosis study, on the task of cell-level binary pathology classification (Fibrotic = Fibroblastic Focus + Fibrosis + Severe Fibrosis, 63.5% prevalence; Control = Minimally Remodeled Alveoli, 36.5%), with inverse-frequency class weights addressing the 1.74:1 imbalance. Of the 343-gene panel, 217 map to a four-level Reactome hierarchy and serve as input features. Cell-cell graphs are constructed per tissue section as a k = 6 nearest-neighbour graph over Euclidean x/y centroid coordinates (1,985,946 edges, zero cross-sample), with edge homophily h = 0.9998 because pathology annotations occupy spatially contiguous polygons.

As for the bulk multi-omics data, we evaluated the pancreatic adenocarcinoma cohort (TCGA-PAAD, N = 175, recurrence rate 46%) released by Zhao et al. as part of MULGONET, on the task of binary recurrence prediction. The cohort provides four omics tracks per patient: mRNA expression, Illumina HumanMethylation450 β-values, and copy-number variation split by sign into amplification (positive log₂ ratios) and deletion (negative log₂ ratios, absolute value) tracks; we used the authors' preprocessed archive directly to isolate model effects from preprocessing variability. The top 1,000 features per modality were selected by a univariate χ² test against the recurrence label and concatenated in the fixed order mRNA, methylation, amplification, deletion to form a 4,000-dimensional patient feature vector — with the caveat, inherited from MULGONET, that selection is performed once on the full cohort and therefore leaks label information into the feature panel; because every model receives the same panel, this does not bias the relative ranking between architectures. The pathway backbone is a Gene Ontology hierarchy restricted to the Biological Process (root GO:0008150) and Molecular Function (root GO:0003674) namespaces treated in parallel, layerized by breadth-first search from each root into five strata. GO terms annotating more than 200 genes were pruned. The retained child-to-parent adjacency at each level defines the sparse connectivity of the corresponding pathway-encoder layer, so the encoder's structure mirrors the ontology exactly, yielding a 29-dimensional pathway embedding (BP + MF concatenated) per patient.

The CITE-seq cohort of 174,753 cells from 41 donors (33 COVID-positive, 8 healthy controls), profiled across two modalities: transcriptome-wide RNA (34,548 genes, SCTransform log-normalized) and a 224-antibody protein panel (ADT, log-normalized), with cells further split into three lineages — NK, T (CD4+CD8), and B (B+plasmablast) — based on coarse manual cell-type annotations. The prediction task is donor-level binary classification of COVID status from cell-level multi-omic profiles. To impose a biological prior on the encoder, gene and antibody features were mapped onto a three-level Reactome pathway hierarchy: RNA gene symbols were matched directly, while antibody panel names were first resolved to their corresponding HGNC gene symbols via a curated synonym table, after which both modalities were restricted to the intersection with Reactome-annotated genes and merged into a single combined feature set (9,661 features) per lineage. Cell-cell graphs were then constructed independently for each lineage as a k-nearest-neighbour graph (k=5) over the dataset’s pre-computed 50 dimensional principal component embedding, yielding one graph per lineage that serves as input to the pathway-masked graph encoders.

(insert figure here of enesemble architecture)

Encoder and Training params

Pathway-masked graph encoders

All three benchmark tasks use the same underlying encoder design: a stack of pathway-pruned linear blocks whose sparse connectivity mirrors a curated gene–pathway hierarchy (Gene Ontology for bulk data, Reactome for spatial and single-cell data), with the propagation operator inside each block swapped between three variants — a pathway-pruned linear layer with no graph structure (ANN), degree-normalized graph convolution (GCN), and learned graph attention (GAT). The three variants are otherwise identical in every other design choice, isolating the effect of the propagation operator (Results).

General formulation. For a hierarchy with $L$ levels, let ${u_0, u_1, \ldots, u_L}$ denote the layer widths, where $u_0$ is the number of input features (genes or gene/protein pairs mapped onto the hierarchy) and $u_\ell$ for $\ell > 0$ is the number of hierarchy terms retained at level $\ell$ after pruning. At each level transition $\ell \in {0, \ldots, L-1}$, a binary adjacency mask $\mathbf{M}\ell \in {0,1}^{u{\ell+1} \times u_\ell}$ encodes parent–child connectivity in the hierarchy:

$$[\mathbf{M}_\ell]_{ij} = \begin{cases} 1 & \text{if term (or gene) } j \text{ at level } \ell \text{ is a child of term } i \text{ at level } \ell{+}1, \ 0 & \text{otherwise.} \end{cases}$$

Each block's weight matrix $\mathbf{W}\ell \in \mathbb{R}^{u{\ell+1} \times u_\ell}$ is constrained to this sparsity pattern by the Hadamard product $\widetilde{\mathbf{W}}\ell = \mathbf{W}\ell \odot \mathbf{M}_\ell$, implemented as a forward pre-hook (torch.nn.utils.prune.custom_from_mask) so that gradients accumulate only on biologically meaningful entries.

Given an input feature matrix $\mathbf{X} \in \mathbb{R}^{N \times u_0}$ and, for graph variants, a symmetric-normalized adjacency $\hat{\mathbf{A}} = \hat{\mathbf{D}}^{-1/2}(\mathbf{A}+\mathbf{I})\hat{\mathbf{D}}^{-1/2}$ built from a $k$-nearest-neighbor similarity graph, the block update at level $\ell$ is

$$\mathbf{H}_{\ell+1} = \operatorname{Dropout}_p\Bigl(\sigma\bigl(\operatorname{Norm}\bigl(\phi_\ell(\mathbf{H}_\ell)\bigr)\bigr)\Bigr), \qquad \mathbf{H}_0 = \mathbf{X},$$

where $\sigma(\cdot) = \operatorname{ReLU}(\cdot)$ and the block operator $\phi_\ell$ is defined per variant as

$$\phi_\ell^{\text{ANN}}(\mathbf{H}) = \mathbf{H},\widetilde{\mathbf{W}}_\ell^\top + \mathbf{1}_N \mathbf{b}_\ell^\top,$$ $$\phi_\ell^{\text{GCN}}(\mathbf{H}) = \hat{\mathbf{A}},\mathbf{H},\widetilde{\mathbf{W}}_\ell^\top + \mathbf{1}_N \mathbf{b}_\ell^\top,$$ $$\phi_\ell^{\text{GAT}}(\mathbf{H}) = \bigl(\boldsymbol{\alpha}_\ell \odot_{\mathcal{E}} \mathbf{H}\bigr),\widetilde{\mathbf{W}}_\ell^\top + \mathbf{1}_N \mathbf{b}_\ell^\top,$$

with attention coefficients, for head $h \in {1, \ldots, H}$,

$$\alpha_{\ell,ij}^{(h)} = \operatorname{softmax}_{j \in \mathcal{N}(i)}\Bigl(\operatorname{LeakyReLU}\bigl(\mathbf{a}_\ell^{(h)\top}\bigl[\widetilde{\mathbf{W}}_\ell^{(h)} \mathbf{h}_i ,\Vert, \widetilde{\mathbf{W}}_\ell^{(h)} \mathbf{h}_j\bigr]\bigr)\Bigr) \quad (\text{LeakyReLU slope } 0.2),$$

where $\mathcal{N}(i)$ is the $k$-nearest-neighbor set of node $i$ and $\Vert$ denotes concatenation. Attention coefficients are recomputed at every forward pass and are not themselves pathway-masked; only the linear projection $\widetilde{\mathbf{W}}_\ell$ is pruned. When $H > 1$, the pathway mask is repeated across heads at the first layer and expanded block-diagonally at subsequent layers, so every head preserves the hierarchy connectivity constraint.

The three datasets instantiate this formulation with a shared branch structure (single branch, dual ontology branches, or per-lineage branches), dataset-specific hierarchy widths, and dataset-specific normalization and head-count choices, summarized in Table 1 and detailed below.

	Bulk (PAAD)	Spatial	CITE-seq
Hierarchy	Gene Ontology (BP + MF)	Reactome	Reactome
Branches	2 (BP, MF)	1	3 (NK, T, B)
Levels $L$	5	4	3
Input width $u_0$	4,000	217	9,661
GAT heads $H$	1	3	1
Block normalization	BatchNorm	BatchNorm (ANN/GCN); LayerNorm (GAT)	BatchNorm
Graph construction	patient-similarity $k$NN	spatial $k{=}6$ $x/y$-centroid $k$NN	per-lineage $k$NN on 50-PC embedding

Bulk data (PAAD). Following Zhang et al.'s MULGONET design, two parallel branches — one following the GO Biological Process (BP) sub-DAG, one following Molecular Function (MF) — each instantiate the general formulation above with $L=5$, $u_0^{(o)} = 4{,}000$, and

$${u_\ell^{(\text{BP})}}_{\ell=0}^{5} = (4000, 2819, 2252, 534, 114, 19), \qquad {u_\ell^{(\text{MF})}}_{\ell=0}^{5} = (4000, 1428, 522, 172, 37, 10).$$

Both branches use BatchNorm and a single attention head ($H{=}1$) for the GAT variant. After the final block, an unpruned linear head projects each branch to its per-root activation vector, $C^{(\text{BP})} = 19$ and $C^{(\text{MF})} = 10$ dimensions respectively. These are concatenated into a 29-dimensional patient embedding $\mathbf{Z} = [\mathbf{Z}^{(\text{BP})} \Vert \mathbf{Z}^{(\text{MF})}]$ and passed through an unpruned, fully dense two-layer fusion MLP (hidden width $d_h$, dropout $p=0.2$) to a single logit, trained end-to-end with binary cross-entropy against the recurrence label:

$$\hat{\mathbf{y}} = \operatorname{sigmoid}\bigl(\operatorname{Dropout}_p(\operatorname{ReLU}(\mathbf{Z}\mathbf{W}_{\text{fus},1}^\top + \mathbf{b}_{\text{fus},1})),\mathbf{W}_{\text{fus},2}^\top + b_{\text{fus},2}\bigr).$$

Spatial transcriptomics. All three variants share a single branch instantiating the general formulation with $L=4$, $u_0 = 217$ Reactome-mapped genes, and ${u_\ell}_{\ell=0}^{4} = (217, 293, 129, 52, 19)$, over a per-sample $k{=}6$ spatial $k$NN graph on Euclidean $x/y$ centroid coordinates ($N=330{,}991$ cells; 1,985,946 edges, zero cross-sample edges). Unlike the other two datasets, the GAT variant here uses $H=3$ heads with concatenation, widening the effective hidden widths to $(217, 879, 387, 156, 57)$ (607,691 parameters, versus 110,593 shared by ANN and GCN), and replaces BatchNorm with LayerNorm to prevent the distributional-shift collapse otherwise caused by BatchNorm's running statistics mismatching test-fold activation distributions under GAT's input-dependent attention weights. After the final block, an unpruned linear head projects directly to $C=2$ class logits, trained full-batch with inverse-frequency-weighted cross-entropy (weights recomputed per fold). Transductive and inductive evaluation differ only in whether $\hat{\mathbf{A}}$ includes validation/test edges during training; five-fold patient-stratified cross-validation is repeated across three seeds with a per-fold-refit StandardScaler.

CITE-seq. Each of the three cell lineages (NK, T, B) is served by an independent branch instantiating the general formulation with $L=3$, $u_0 = 9{,}661$ Reactome-mapped genes/proteins, and ${u_\ell}{\ell=0}^{3} = (9{,}661, 400, 113, 24)$ (mask densities $0.61%$, $0.89%$, $4.20%$ across the three level transitions), over a per-lineage $k$NN graph on the dataset's 50-dimensional PC embedding. The GAT variant uses a single head ($H{=}1$, unlike the spatial encoder's $H{=}3$). After the final block, an unpruned linear head projects each lineage's root-pathway activations to a $d=32$-dimensional per-cell embedding, $\mathbf{Z}^{(m)} = \mathbf{H}L^{(m)}\mathbf{W}{\text{emb}}^{(m)\top} + \mathbf{b}{\text{emb}}^{(m)\top}$.

Because prediction is at the donor level while the encoder operates on cells, an attention-pooling step aggregates each lineage's cell embeddings per donor. For donor $d$, cell set $\mathcal{C}_d^{(m)} = {i : \text{donor}(i)=d}$ in lineage $m$:

$$\mathbf{u}_i^{(m)} = \tanh\bigl(\mathbf{W}_{\text{pool}}^{(m)} \mathbf{z}_i^{(m)} + \mathbf{b}_{\text{pool}}^{(m)}\bigr), \qquad w_i^{(m)} = \frac{\exp(\mathbf{v}^{(m)\top}\mathbf{u}_i^{(m)})}{\sum_{j \in \mathcal{C}_d^{(m)}} \exp(\mathbf{v}^{(m)\top}\mathbf{u}_j^{(m)})}, \qquad \mathbf{e}_d^{(m)} = \sum_{i \in \mathcal{C}_d^{(m)}} w_i^{(m)}\mathbf{z}_i^{(m)}.$$

Donors lacking cells in a given lineage receive a zero embedding and a mask bit $m_d^{(m)}=0$. The three pooled lineage embeddings and their mask bits are concatenated into a 99-dimensional donor representation $\mathbf{r}_d = [\mathbf{e}_d^{(\text{NK})}\Vert\mathbf{e}_d^{(\text{T})}\Vert\mathbf{e}_d^{(\text{B})}\Vert m_d^{(\text{NK})}\Vert m_d^{(\text{T})}\Vert m_d^{(\text{B})}]$ and passed through an unpruned, fully dense two-layer fusion MLP ($d_h=128$, dropout $p_f=0.3$) to $C=2$ class logits, trained with softmax cross-entropy against donor COVID status.

END

Stats:

Spatial

Edge Homophily: 0.9998

Node Homophily: 0.9998

Bulk (Mulgo)

Edge: 0.5907

Node: 0.6116

With cosine z-score graph: see improvement in performance

h_edge: 0.641

h_node : 0.279

Citeseq

# nk edge homophily: 0.6941, node homophily: 0.6979

# t edge homophily: 0.7696, node homophily: 0.7697

# b edge homophily: 0.8269, node homophily: 0.8217

Plots:

1x Table for reuslts. 3x5 . rows: datasets. Columns: AUROC

1x Ensemble workflow/diagram

3x violin plots for the three models across three data types. Do thresholded F1 as well ?

3x parallel Precision-recall curve.

Optional: . The neighborhood class-distribution heatmap for the three datasets.

TODO:

• Write Results paragraphs on performance. Generate plots where necessary. Figure out what the plots should be

Goal: Have a draft of the paper finished by the end of the week.

Technical TODO:

• Clean up all code not related to library

• Add GO integration (Do)

• Allow for protein synonyms to be passed in model (Do next)

• Add visual plotting,

• Add SHAP?

• Option to choose final hyperparameter (Do)

• Pruning pathways for which we don’t have feature cells. Must be added

• Allow for regression, clean the library up

• Allow for Ann Data to graph data to setup its own structure

• Modify the ensemble class to take the encoders as input. Fucking hate the api right now.

Interpretability features would be:

1. SHAP

2. Integrated gradients

3. Sankey diagram

Experimental TODO:

1. run Covid Sequence experiment. Record results. Save Model and evaluation results including plots. Write code for training and testing.

2. Do exact samel experiment with Mulgo

3. Figure out the new data and run on that as well.

Additional experiments:

Use integrated gradients and

Writing TODO:

1. Keep track of citations

2. Write paper :)

Hypothesis:

• Bio-informed GAT and GCN are better than ANN:
  Reasoning:

Citeseq cleanup:

There’s a 10 percent difference between BINN GCN and BINN ANN’s accuracy.
Need to also run on BINN GAT to record in the paper. Put a seed with that in the paper.

Comparison against MULGONET-mode baseline

Model Clean library (this session) MULGONET-mode (prior session) Δ ANN 0.6525 ± 0.0777 0.8026 ± 0.0798 −0.150 GCN 0.6806 ± 0.0843 0.8815 ± 0.0274 −0.201 GAT 0.7196 ± 0.0612 0.8595 ± 0.0401 −0.140

Will include both methods in the paper.

Data + Results

Single cell:

Used Reactome hierarchy
Has an attention encoder to begin in the loop.
CITE-seq data from COVID PBMCs (Peripheral Blood Mononuclear Cells). -> Round nucleas cells in blood.
Different genes get expressed under differing scenarios like covid / no covid.

Results:

Bulk data

Mulgonet: https://www.sciencedirect.com/science/article/pii/S2667325825000147

Here are the exact measurements from the paper-identical run (PAAD, paper’s exact GitHub data, 175 samples, 4 omics modalities, 5-fold stratified CV, 150 epochs):

Model AUPR AUROC F1 Accuracy Precision Recall ANN 0.8026 ± 0.0798 0.8042 ± 0.0838 0.6599 ± 0.1023 0.6571 ± 0.0404 0.6201 0.7831 GCN 0.8815 ± 0.0274 0.8771 ± 0.0142 0.7373 ± 0.0656 0.7657 ± 0.0457 0.7625 0.7265 GAT 0.8595 ± 0.0401 0.8581 ± 0.0268 0.7283 ± 0.0402 0.7600 ± 0.0428 0.7767 0.6912