Journal Title Here DOI HERE Advance Access Publication Date: Day Month Year
Paper

Bontonou et al.

[$\ast$]Corresponding authors<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
0Year 0Year 0Year

# A Comparative Analysis of Gene Expression Profiling by Statistical and
Machine Learning Approaches

Myriam Bontonou0000-0002-0010-5457  Anaïs Haget  Maria Boulougouri  Benjamin
Audit0000-0003-2683-9990  Pierre Borgnat0000-0003-4536-8354  Jean-Michel
Arbona0000-0001-6166-9056  Laboratoire de Biologie et Modélisation de la
Cellule, ENS de Lyon, CNRS, UMR5239, Inserm U1293, Univ Lyon, France  Univ
Lyon, ENS de Lyon, CNRS, Laboratoire de Physique, F-69342 Lyon, France  LTS2
laboratory, EPFL, Lausanne, Switzerland

(2024; 2024; Date; Date; Date)

###### Abstract

Many machine learning models have been proposed to classify phenotypes from
gene expression data. In addition to their good performance, these models can
potentially provide some understanding of phenotypes by extracting
explanations for their decisions. These explanations often take the form of a
list of genes ranked in order of importance for the predictions, the highest-
ranked genes being interpreted as linked to the phenotype. We discuss the
biological and the methodological limitations of such explanations.
Experiments are performed on several datasets gathering cancer and healthy
tissue samples from the TCGA, GTEx and TARGET databases. A collection of
machine learning models including logistic regression, multilayer perceptron,
and graph neural network are trained to classify samples according to their
cancer type. Gene rankings are obtained from explainability methods adapted to
these models, and compared to the ones from classical statistical feature
selection methods such as mutual information, DESeq2, and EdgeR.
Interestingly, on simple tasks, we observe that the information learned by
black-box neural networks is related to the notion of differential expression.
In all cases, a small set containing the best-ranked genes is sufficient to
achieve a good classification. However, these genes differ significantly
between the methods and similar classification performance can be achieved
with numerous lower ranked genes. In conclusion, although these methods enable
the identification of biomarkers characteristic of certain pathologies, our
results question the completeness of the selected gene sets and thus of
explainability by the identification of the underlying biological processes.

###### keywords:

Machine learning, Explainability, Gene expression data, Logistic regression,
Multilayer perceptron, Graph neural network, Differential expression analysis,
Over-representation analysis

##  1 Introduction

A phenotype results from a complex cascade of molecular processes possibly
involving hundreds or thousands of genes. Understanding the connections
between phenotype and gene expression is crucial to better identify mechanisms
implicated in diseases such as cancer.

Statistical methods like EdgeR robinson2010edger and DESeq2 love2014moderated
, have been pivotal in analysing gene expression levels. They identify genes
that are expressed differentially between phenotypes and suggest to select
them as markers of these phenotypes. However, recent Machine Learning (ML)
methods, having proved their usefulness in other fields Goodfellow-et-al-2016
, open a promising avenue to enhance diagnostic accuracy and unravel complex
biological processes. This has been showcased for identifying cancer subtypes
golub1999molecular ; hoadley2018cell ; parker2009supervised and for training
diagnostic tools on large datasets leng2022benchmark ; hanczar2022assessment ,
leading to the identification of discriminating sets of genes
li2017comprehensive ; jacquet2023aberrant . Various models, such as
$k$-nearest neighbours li2017comprehensive , support vector machines
rohimat2022implementation , deep neural networks ahn2018deep and graph neural
networks ramirez2020classification , have been successfully trained to
classify samples from their gene expression profiles. Discriminating gene sets
can be determined either before, during, or after training
mahendran2020machine ; alharbi2023machine . Integrating these ML methods in
the analysis of gene expression levels suggests a paradigm shift in these
analyses: interpretability (or explainability) techniques for ML molnar2022
offer ways to explain the decisions of the models, and provide scores that
allow to rank and possibly select relevant genes to explain specific
phenotypes. However, there is a significant gap between the selection of
discriminating gene sets and our understanding of the biological processes
involved in the development of certain phenotypes, because distinct gene sets
can yield similar classification performance li2017comprehensive .

Here, we study whether ML models and associated scores of explainability can
unveil novel, biologically relevant, molecular signatures. The signatures
obtained from several statistical and ML methods are studied in terms of their
classification performance and their biological relevance. Biological
relevance is evaluated through over-representation analysis (ORA)
khatri2012ten , which computes the overlap of the top-ranked genes with
established gene sets representing various biological processes. In this
context, gene expression profiling refers to the creation of a profile of
genes ranked according to their relevance for understanding (and possibly
predicting) phenotypes. To elucidate the specific insights that ML can offer,
the top ranked genes selected by ML methods are compared with those selected
by EdgeR robinson2010edger and DESeq2 love2014moderated methods for
differential gene expression. Several classifiers, including logistic
regressions, shallow neural networks leveraging more complex relationships,
and graph-based neural networks exploiting gene interactions, are trained to
differentiate cancer tissue samples. Genes rankings are derived from the ML
classifiers using the integrated gradients explainability method (IG)
sundararajan2017axiomatic .

Experiments are carried out on existing gene expression datasets from The
Cancer Genome Atlas (TCGA)111www.cancer.gov/tcga, the Therapeutically
Applicable Research to Generate Effective Treatments programme
(TARGET)222www.cancer.gov/ccg/research/genome-sequencing/target and the
Genotype-Tissue Expression project (GTEx) lonsdale2013genotype .

The results, in line with existing literature, identify small gene sets
maintaining classification performance. Still, the question of the relevance
of genes selected through ML models explainability remains unsolved. We
observe that the top-ranked genes identified by ML methods vary, and differ
significantly from those identified by statistical methods – presenting
themselves significant variations. Interestingly, a classifier trained on
genes selected by DESeq2 or mutual information outperforms one trained on
genes specifically selected for it. Over-representation analysis reveals
diverse biological processes, sometimes specific to a single method or a
paired methods, suggesting different facets of explainability. In simple
cases, the ranks obtained by explaining ML models predictions correlate well
with t-statistic ranks comparing gene means of two classes, which is aligned
with our intuitive understanding of the problem.

In summary, we undertake a comprehensive analysis of gene expression profiling
using emerging ML techniques by questioning classification efficiency and
biological relevance. Comparisons with traditional statistical methods provide
insights into the evolving landscape of molecular signature identification.
The code and data are accessible at https://github.com/mbonto/XAI_in_genomics.

##  2 Materials and methods

###  2.1 Setting

Let’s consider a gene expression dataset containing $N$ training data samples
of class $c\in\llbracket 1;C\rrbracket$. The class reflects the phenotypic
state of a tissue. Here, it is a cancer type, a cancer subtype, or a healthy
type. A data sample is a feature vector $\mathbf{{x}}\in\mathbb{R}^{G}$
containing the average expression of $G$ genes within a tissue sample. The
true class $c$ of a sample is represented by a one-hot vector
$\hat{\mathbf{{y}}}\in\mathbb{R}^{C}$ having value 0 everywhere except for
$\hat{\mathbf{{y}}}_{c}=1$. The Euclidean norm is written
$\left\lVert\cdot\right\rVert$.

###  2.2 Machine learning models for classification

To solve a classification task, a supervised model
$f:\mathbb{R}^{G}\mapsto\mathbb{R}^{C}$ learns to map the features
$\mathbf{{x}}$ of a data sample to a vector $\mathbf{{y}}\in\mathbb{R}^{C}$
whose coefficients $\mathbf{{y}}_{c}$ represent the probability of belonging
to class $c$. Three types of models are used in the present work: logistic
regression (LR), multilayer perceptron (MLP), and graph neural network (GNN).

####  2.2.1 Logistic regression (LR)

LR models with L2 or L1 regularisation for both binary ($C=2$ classes) and
multi-class ($C>2$) classification problems are considered. In the binary
case, a sigmoid function is applied to compute class probabilities
bishop2006pattern . Given the parameters $\mathbf{{W}}\in\mathbb{R}^{1\times
G}$, $b\in\mathbb{R}$ and $\sigma(z)=1/(1+\exp(\scalebox{0.75}[1.0]{$-$}z))$
the sigmoid function, the model is
$\mathbf{{y}}=\sigma(\mathbf{{W}}\mathbf{{x}}+b)$. In the multi-class case, a
softmax function is used instead of $\sigma$ for the probabilities
bishop2006pattern . Given the parameters $\mathbf{{W}}\in\mathbb{R}^{C\times
G}$ and $\mathbf{{b}}\in\mathbb{R}^{C}$,
$\mathbf{{y}}=\text{softmax}(\mathbf{{W}}\mathbf{{x}}+\mathbf{{b}})\text{ with
softmax}(\mathbf{{z}}_{g})=\frac{\exp(\mathbf{{z}}_{g})}{\sum_{f}\exp(\mathbf{{z}}_{f})}\,.$
Two versions of LR, LR+L2 and LR+L1, are considered. LR+L2 is trained with a
L2 penalty, a regularisation term equal to the squared values of the
parameters. LR+L1 uses a L1 penalty term equal to the absolute values of the
parameters; it has the advantage of sparsity. As the number of non-zero
parameters is minimised, the interpretability of potential candidate genes
relevant for phenotypes is enhanced. These regularisation techniques both
prevent the model from over-fitting to the training data, and improve the
generalisation performance bishop2006pattern .

####  2.2.2 Multilayer perceptron (MLP)

A MLP is a neural network containing a series of fully connected layers
followed by a LR Goodfellow-et-al-2016 . Formally, $G^{[l]}$ is the number of
hidden features after layer $l$ with $G^{[0]}=G$. In the binary case,
$f(\mathbf{{x}})=\sigma(\mathbf{{b}}^{[L]}+\mathbf{{W}}^{[L]}\text{\small{ReLU}}(\dots\text{\small{ReLU}}(\mathbf{{b}}^{[1]}+\mathbf{{W}}^{[1]}\mathbf{{x}})))$.
The parameters learned within layer $l$ are
$\mathbf{{W}}^{[l]}\in\mathbb{R}^{G^{[l]}\times G^{[l-1]}}$ and
$\mathbf{{b}}^{[l]}\in\mathbb{R}^{G^{[l]}}$. ReLU is the Rectified Linear Unit
function. For the multi-class case, the sigmoid of the last layer is replaced
by softmax. Here, the MLP architecture also includes a batch normalisation
function to each layer, which stabilises training and improves generalisation
Goodfellow-et-al-2016 . In the experiments, shallow MLP with 1 or 2 layers are
used (see hyperparameter selection).

####  2.2.3 Graph neural network (GNN)

A GNN zhou2020graph incorporates gene pairs relationships using a graph
structure. Following ramirez2020classification , a graph is created by
connecting co-expressed genes through Pearson correlations. This generates a
graph $\mathcal{G}$, with nodes $\mathcal{V}=\left\\{1,\dots,G\right\\}$
representing genes and edges $\mathcal{E}$ denoting relationships between gene
pairs. The edge weights (i.e., thresholded correlations) are stored in the
adjacency matrix $\mathbf{{A}}\in\mathbb{R}^{G\times G}$ of $\mathcal{G}$. The
applied threshold is discussed in the hyperparameter selection subsection. As
in ramirez2020classification , the GNN architecture includes 1 or 2 neural
network layers followed by a LR. Each neural network layer comprises a graph
convolutional layer DBLP:conf/iclr/KipfW17 followed by a graph coarsening
layer based on the Graclus algorithm dhillon2007weighted .

Formally, $F[l]$ represents the number of features associated with a gene
after layer $l$. In the initial layer, the gene expression value is the only
feature ($F[0]=1$). The first graph convolutional layer is
$\text{ReLU}(\mathbf{{b}}^{[1]}+\tilde{\mathbf{{D}}}^{-\frac{1}{2}}\tilde{\mathbf{{A}}}\tilde{\mathbf{{D}}}^{-\frac{1}{2}}\mathbf{{x}}\mathbf{{W}}^{[1]})$
with the data sample $\mathbf{{x}}\in\mathbb{R}^{G\times 1}$, the normalised
adjacency matrix $\tilde{\mathbf{{A}}}=\mathbf{{A}}+\mathbf{{I}}$ and
$\tilde{\mathbf{{D}}}$ the degree matrix of $\tilde{\mathbf{{A}}}$
DBLP:conf/iclr/KipfW17 . The same parameters
$\mathbf{{W}}^{[1]}\in\mathbb{R}^{1\times F[1]}$ and
$\mathbf{{b}}\in\mathbb{R}^{F[1]}$ are applied to all genes. The graph
coarsening layer dhillon2007weighted combines pairs of adjacent nodes into a
single node, preserving maximal features (max pooling). New edges result from
the union of previous edges, with associated weights summed.

###  2.3 Gene selection and explainability for gene profiling

The objective is to identify biologically significant molecular signatures
using gene expression data. Several methods assign a score
$\boldsymbol{\phi}_{g}$ to each gene $g$ (its expression level being an input
feature) indicating its relative importance for the sample phenotype class.
These methods are categorised into three groups. (i) Filter methods rank genes
based on a statistics measuring the amount of information they contain;
filtering is done independently of any specific classifier. (ii) Embedded
methods rank genes based on a score derived from a classifier; gene ranking is
integrated into the training process of the classifier. (iii) Post-hoc methods
rank genes based on a score computed after the training of a classifier; gene
ranking occurs as a separate step after the classifier is trained.

####  2.3.1 Filter method - Variance (VAR)

Genes can be ranked in decreasing order of their variance across samples in
the dataset; filtering out genes with low variance is often standard practice
when preprocessing a dataset. By focusing on the gene expression data only and
ignoring their corresponding phenotypes, this ranking method disregards true
classes. Let $\overline{\mathbf{{x}}}\in\mathbb{R}^{G}$ be the average on all
the training samples. Then, the score is:

$\boldsymbol{\phi}_{g}^{\text{VAR}}=\frac{1}{N}\sum_{n=1}^{N}(\mathbf{{x}}_{g}^{n}-\overline{\mathbf{{x}}}_{g})^{2}\,.$
(1)

####  2.3.2 Filter method - Principal component analysis (PCA)

Genes can be ranked according to their contribution to the main directions of
data variability. The higher the magnitude of a gene’s coefficient in the
principal components, the more significant it is jolliffe2002principal . PCA
focuses on the variability of the dataset without considering the classes. As
large proportion of this variability is caught by the first principal
component $\mathbf{{v}}^{1}$, i.e. associated with the highest eigenvalue of
the covariance $\mathbf{{\Sigma}}\in\mathbb{R}^{G\times G}$, only this one is
kept here::

$\boldsymbol{\phi}_{g}^{\text{PCA}}=\lvert\mathbf{{v}}^{1}_{g}\rvert\,.$ (2)

####  2.3.3 Filter method - Mutual information (MI)

To account for higher order terms in probabilities, genes are ranked in
decreasing order of the mutual information $I(X_{g};Y)$ shared between the
expression $X_{g}$ of a gene $g$, and the classes represented by a discrete
variable $Y$ ross2014mutual . Then:

$\boldsymbol{\phi}_{g}^{\text{MI}}=I(X_{g};Y)\,.$ (3)

####  2.3.4 Filter method - Differential expression (EdgeR, DESeq2)

Popular bioinformatics tools are used to identify differentially expressed
genes between experimental conditions costa2017rna . Once the expression
distribution of a gene is modelled, a test for differential expression is
performed and yields a p-value adjusted to account for multiple testing.

In this study, we consider EdgeR robinson2010edger and DESeq2
love2014moderated , two methods of this category that rely on different
testing procedures. Genes are then ranked according to their adjusted
p-values; low adjusted p-values indicating a high statistical significance for
differential expression. When more than two classes are considered, a test is
performed for each possible pair of classes. The genes are ranked based on the
minimal p-value obtained across all these tests. This ranking strategy
reflects the presence of a significant difference for at least one pair of
classes. The scores are:

$\boldsymbol{\phi}_{g}^{\text{EdgeR, DESeq2}}=-\text{log}_{10}(\text{adjusted
p-value}^{\text{EdgeR, DESeq2}}_{g})\,.$ (4)

####  2.3.5 Embedded method - Magnitude of the LR weights

For LR, the magnitude of the weight $\mathbf{{W}}_{g}$ associated to a gene
$g$ reflects its impact on the classification, as a larger one signifies a
more substantial influence of the associated gene. If it has positive value,
it indicates that over-expression favours the class, while negative value
suggests that under-expression does so. For binary classification, an
intrinsic score is directly the magnitude of this weight molnar2022 . As gene
expressions are standardised, this metric ensures comparability of magnitudes
across genes. The score for binary LR is:

$\boldsymbol{\phi}_{g}^{\text{LR (weight)}}=|\mathbf{{W}}_{g}|\,.$ (5)

For multi-class, several linear functions (one for each class) are
simultaneously learned. By analogy, genes are ranked based on the average of
the absolute values of the parameters associated with a gene, hence:

$\boldsymbol{\phi}_{g}^{\text{LR
(weight)}}=\frac{1}{C}\sum_{c=1}^{C}\lvert\mathbf{{W}}_{cg}\rvert\,.$ (6)

####  2.3.6 Post-hoc method - Integrated Gradients (IG)

Neural networks are not directly interpretable because of their functional
complexity molnar2022 . To highlight the individual features that impact the
most the decision of a neural network on a particular example, several
explainability methods have been developed lundberg2017unified . Here, the
integrated gradients method (IG) sundararajan2017axiomatic is chosen to
explain the decisions made by various gradient-based ML models, specifically
LR, MLP, and GNN. IG is a gradient-based technique that is used a lot in the
literature, in particular for its computational efficiency.

IG assigns a score function
$\boldsymbol{\phi}_{g}^{\text{local}}(\mathbf{{x}})$ to each gene $g$ in each
sample $\mathbf{{x}}$ to represent its importance for the model’s decision.
The computation of the scores
$\boldsymbol{\phi}_{g}^{\text{local}}(\mathbf{{x}})$ contrasts the prediction
for a sample $\mathbf{{x}}$ with a reference sample $\mathbf{{x}}^{\prime}$
(baseline) according to:

$\boldsymbol{\phi}_{g}^{\text{local}}(\mathbf{{x}})=(\mathbf{{x}}_{g}-\mathbf{{x}}_{g}^{\prime})\int_{\alpha=0}^{1}\frac{\partial
f_{c}(\mathbf{{z}})}{\partial\mathbf{{x}}_{g}}\bigg{|}_{\mathbf{{z}}=\mathbf{{x}}^{\prime}+\alpha(\mathbf{{x}}-\mathbf{{x}}^{\prime})}d\alpha\,.$

Here, $f_{c}(\mathbf{{z}})$ is the output of the model for class $c$ at input
$\mathbf{{z}}$. To aggregate the decisions of a model, a global score
$\boldsymbol{\phi}_{g}^{\text{class}}$ for the importance of a gene $g$ within
a class $c$ is derived using training samples that are correctly classified.
If $N_{\text{class}}$ is the number of data samples of class $c$ used to train
a model, then:

$\boldsymbol{\phi}_{g}^{\text{class}}=\frac{1}{N_{\text{class}}}\sum_{n=1}^{N_{\text{class}}}\frac{\lvert\boldsymbol{\phi}_{g}^{\text{local}}(\mathbf{{x}}^{n})\rvert}{\left\lVert[\boldsymbol{\phi}_{1}^{\text{local}}(\mathbf{{x}}^{n}),\dots,\boldsymbol{\phi}_{G}^{\text{local}}(\mathbf{{x}}^{n})]\right\rVert}$
(7)

The normalisation ensures comparability across different samples. When
multiple classes are studied, a global metric is derived by averaging the
importance scores:

$\boldsymbol{\phi}_{g}^{\text{IG}}=\frac{1}{C}\sum_{c}\boldsymbol{\phi}_{g}^{\text{class}}\,.$
(8)

The studied classes and the baselines used for each classification task
considered in this article are detailed in the feature selection methods
subsection of the Materials Section.

####  2.3.7 Prediction Gaps for IG

One can obtain more insights into the significance of the IG ranking in the
classification process by studying the impact of progressively masking genes
on the model’s predictions. This masking, referred to as experiment 0 in the
following is performed by replacing the gene expression values in the original
sample $\mathbf{{x}}$ with the corresponding values from a reference sample
$\mathbf{{x}}^{\prime}$. For a sample $\tilde{\mathbf{{x}}}_{m}$ with $m$
masked variables, the prediction gap (PG) agarwal2022openxai is:

$\text{PG}=\frac{1}{G}\sum_{m=1}^{G}\frac{\max(f_{c}(\mathbf{{x}})-f_{c}(\tilde{\mathbf{{x}}}_{m}),0)}{f_{c}(\mathbf{{x}})}\;.$

Genes can be masked using different orders on the predefined rankings
$\boldsymbol{\phi}$: in descending order of importance, the metric is called
PGI (I for important masked first); in ascending order of importance, it is
called PGU (U for unimportant masked first). Fig. 1 shows a sketch of what to
expect of these metrics. Note that the transition between the original
prediction of the data sample and the prediction of the baseline can be less
smooth. If the rankings are good, PGU estimates the fraction of irrelevant
features, while PGI calculates the fraction of important features that are not
necessary for the models.

###  2.4 Assessing the genes selection methods

As the different selection methods generate distinct gene expression profiles,
it is interesting to compare them. Rankings from classifier-based methods
derive from learning, whose stability depends on the optimisation process. To
account for variability, the ML models (LR, MLP, GNN) are trained 10 times
with different initialisations, which generates 10 ranking replicates. In
contrast, VAR, PCA, MI, EdgeR, and DESeq2 each generate a unique ranking per
dataset.

####  2.4.1 Comparing the top-ranked genes across methods

For each dataset, heatmaps depicting the similarity of gene sets obtained by
different methods are displayed (e.g., in Fig. 2). The lower / upper
triangular part of the heatmap displays the percentages of common top 10 / 100
genes across different methods. When there are several ranking replicates, the
percentages are averaged over all possible pairs. The diagonal of the heatmap
shows the average percentage of common genes among the top 100 genes coming
from ranking replicates of a method. For methods with a unique ranking, this
percentage is 100%, by design.

Dataset (source) | Task | # classes | # samples (min/max per class) | # genes
---|---|---|---|---
PanCan (legacy TCGA) | Tumour types | 33 | 9853 (36/1095) | 15401
BRCA (GDC TCGA) | Healthy vs Tumour | 2 | 1210 (113/1097) | 13946
BRCA-pam (legacy TCGA) | PAM50 classes | 5 | 916 (67/421) | 13896
ttg-breast (TCGA TARGET GTEx) | Healthy vs Tumour | 2 | 1384 (292/1092) | 14373
ttg-all (TCGA TARGET GTEx) | Healthy vs Tumour | 2 | 17600 (8130/9470) | 14368
Table 1: Description of the various datasets and their associated
classification tasks.

####  2.4.2 Classification performance investigation

The analysis investigates the minimum number of genes needed for a model to
achieve performance similar to using all genes. Two types of experiments are
conducted. In experiment 1, each model is trained using only the $n_{H}$
highest-ranked genes, with $n_{H}$ from 1 to 1000. In experiment 2, each model
is trained using the $n_{L}$ lowest-ranked genes, $n_{L}$ from 1 to 1000. For
a model that produces several ranking replicates, the performance is averaged
across all replicates. These experiments aim to determine whether (1) the
$n_{H}$ top-ranked genes contain enough information to discriminate the
classes, and (2) whether the $n_{L}$ lowest-ranked genes are similarly
informative.

###  2.5 Biological function investigation

After identifying the 100 most important genes for each method, an over-
representation analysis (ORA) khatri2012ten is performed using gene sets from
the Molecular Signatures Database Human collections liberzon2011molecular
333www.gsea-msigdb.org/gsea/msigdb/human/annotate.jsp. These gene sets,
derived from the analysis of large repositories of gene expression data,
encompass a variety of biological processes and pathological conditions.
Specifically, the overlap of each gene list is calculated with the H hallmark
gene sets liberzon2015molecular , the C2 curated gene sets (including
canonical pathways), the C4 computational gene sets (curated from cancer
expression data), the C5 ontology gene sets, the C6 oncogenic signature gene
sets and the C7 immunologic signature gene sets. Results are filtered
according to significance ($\text{false discovery rate q-value}\leq 0.05$).
For each method, the top 10 over-represented gene sets are displayed. If any
of these gene sets appears among the top 100 significant gene sets identified
by another method, it is flagged. For ML-based methods, only one ranking is
considered.

Most important variables masked firstLess important variables masked
first$0$$10$$20$$30$$40$$50$$60$$70$$80$$90$$100$$0$$50$$100$PGUPGI$\frac{f_{c}(\mathbf{{x}})-f_{c}(\mathbf{{x}}^{\prime})}{f_{c}(\mathbf{{x}})}$Proportion
of masked variables (%)Prediction error(%) Figure 1: Illustration of PGI and
PGU: prediction gaps on important (PGI) and unimportant (PGU) features, for an
example $\mathbf{{x}}$ in a class $c$.

###  2.6 Materials for the experiments

####  2.6.1 Datasets

Bulk RNA-sequencing (RNA-seq) is a powerful technique for quantifying gene
expression levels in tissue samples. The datasets analysed in this study
originate from the TCGA Research Network, the TARGET initiative and the GTEx
project lonsdale2013genotype . They are accessible through the Genomic Data
Commons grossman2016toward 444portal.gdc.cancer.gov or the Xena browser
goldman2020visualizing . TCGA tomczak2015review encompasses gene expression
data from various human tumours and surrounding normal tissues, available in
legacy and Genomic Data Commons (GDC) versions. Notable distinctions between
the two versions include the use of different human reference genomes (hg19
for legacy, hg38 for GDC) gao2019before . The Pan-Cancer Atlas, derived from
legacy TCGA, explores 33 tumour types hoadley2018cell . The TCGA TARGET GTEx
dataset (TTG) integrates gene expression data from adult cancer tissues
(TCGA), pediatric cancer tissues (TARGET), and non-diseased tissues (GTEx),
all processed using a unified pipeline vivian2017toil .

This article focuses on classical classification problems, such as PanCan
ramirez2020classification or BRCA-pam rhee2018hybrid ; choi2023mobrca . The
datasets, detailed in Table 1, involve a limited number of samples for a large
number of features (gene expression). Gene expression is measured by counting
reads (fragment of RNA molecules) in a data sample. Preprocessing steps
include the removal of genes with NaN values, those exhibiting a maximal
expression level of 0 across the dataset, and low-expressed genes with fewer
than 5 counts in over 75% of samples for each class. In the case of ttg-all,
27 samples are excluded due to null expression levels in more than 75% of the
genes. Normalisation involves scaling total counts in each sample to $10^{6}$,
with subsequent $\text{log}_{2}$ transformation of the normalised counts. All
datasets are available for download from our GitHub repository.

####  2.6.2 Hyperparameter selection for ML methods

The data is split randomly with 60% of the samples used for training models
and 40% for testing. The hyperparameters of the models are selected with a
grid search using a 4-split cross-validation on the training data. The
obtained parameters are in Table S1. For the graph, the correlation threshold
is set such that the number of edges is $k\times G$, hence controlling its
density.

Dataset | LR+L1 | LR+L2 | MLP | GNN
---|---|---|---|---
PanCan | $95.0$ | $94.3$ | $94.3\pm 0.3$ | $92.1\pm 0.4$
BRCA | $99.7$ | $98.5$ | $99.5\pm 0.4$ | $98.9\pm 0.6$
BRCA-pam | $92.3$ | $90.7\pm 0.2$ | $87.4\pm 1.8$ | $87.1\pm 1.4$
ttg-breast | $99.7$ | $99.2$ | $99.4\pm 0.3$ | $99.1\pm 0.1$
ttg-all | $99.5$ | $99.5$ | $99.6$ | $99.4\pm 0.1$
Table 2: Classification performance measured by balanced accuracy (%).
Standard deviations are computed from 10 replicates. They are not reported
when below 0.05.

####  2.6.3 Training process for ML methods

Parameters of the models are learned through gradient descent, by minimising
the cross-entropy loss function (plus regularisation). LR+L2 and LR+L1 are
trained with the SAGA solver from scikit-learn with a maximum of 1000
iterations. MLP is trained for 25 epochs using PyTorch’s SGD optimiser
paszke2019pytorch with a momentum of 0.9 and a weight decay of 0.0001. The
initial learning rate is 0.1, decreasing to 0.01 after 13 epochs and to 0.001
after 23 epochs. GNN is trained for 15 epochs using PyTorch’s Adam optimiser
DBLP:journals/corr/KingmaB14 with a weight decay of 0.0001. The initial
learning rate of 0.01 decreases to 0.001 after 8 epochs and to 0.0001 after 14
epochs. To address training variability, the learning process for each model
is repeated ten times, using distinct initialisation seeds. This results in
ten different rankings per model.

####  2.6.4 Feature selection methods

Gene rankings are generated using training samples only. All methods are
implemented in Python except DESeq2 and EdgeR which are based on R packages.
VAR, PCA and MI are implemented using scikit-learn scikit-learn . The code to
include DESeq2 and EdgeR into a python script is inspired from
clarke2021appyters . IG is implemented using the torch package called captum
kokhlikyan2020captum . The IG scores are computed on the tumour class for
BRCA, ttg-breast and ttg-all and on the tumour subtype classes for BRCA-pam.
In these cases, the baseline $\mathbf{{x}}^{\prime}$ with respect to which the
scores are computed is the average of the normal training samples. The scores
are computed on all classes for PanCan with respect to the average of the
training samples.

##  3 Results

###  3.1 Gene expression data is informative on phenotypes

Dataset | LR (non-zero for LR+L1) | MLP | GNN
---|---|---|---
PanCan | 508266 (13575) | 308773 | 969075
BRCA | 13947 (99) | 279001 | 26187
BRCA-pam | 69485 (390) | 278085 | 118484
ttg-breast | 14374 (206) | 287541 | 12501
ttg-all | 14369 (7376) | 576601 | 27456
Table 3: Numbers of parameters learned by each model (averaged for selected
non-zero parameters for LR+L1, and for GNN).

Four ML models, LR+L1, LR+L2, MLP and GNN, are trained to classify gene
expression data from tissues across several cancer types (PanCan), breast
cancer and healthy surrounding tissues (BRCA), breast cancer and healthy
tissues (ttg-breast), several cancer types and healthy tissues (ttg-all) and
various subtypes of breast cancer (BRCA-pam) (Table 1). The classification
performance is evaluated on data samples that have not been seen during
training. Each model is trained 10 times with a different random
initialisation. Table 2 reports the average balanced accuracy scores,
correcting for class imbalance (average of recall obtained on each class).
Unbalanced accuracies are similar, see Table S2. The numbers of parameters for
the models are in Table 3.

Across all datasets and models, phenotypes are predicted with balanced
accuracy consistently exceeding 95% (except for BRCA-pam). High accuracies
exceeding 99% are even achieved for datasets classifying cancer tissues
against normal tissues. For cancer types and subtypes, the best accuracies are
respectively 95% and 92%. LR+L1 often outperforms other models (at the price
of larger training time, see Table S3), with LR+L2 and MLP following closely.

Good performance suggests also that the models have learned informative
patterns linking gene expression data to phenotypes. As LR+L1 performed
similarly to MLPs and GNNs, it suggests that precise classification can be
achieved without necessarily learning interactions between genes; GNNs do not
appear to benefit much from the graph of correlations.

###  3.2 Top-ranked genes obtained by explainability and statistical gene
selection methods differ significantly

((a)) ttg-breast

((b)) BRCA-pam

Figure 2: Heatmaps showing the percentage of common genes among the top 10
(lower) and top 100 (upper + diagonal) genes selected by each method. More
details are in Materials and methods.

For each dataset, genes are ranked by the methods associated with the ML
models (LR, MLP, GNN) or the statistical filtering methods (VAR, PCA, MI,
EdgeR, DESeq2). Heatmaps are constructed to visualise the similarity of the
top-ranked genes (Fig. 2 for datasets ttg-breast and BRCA-pam and Fig. S1 for
the other datasets).

For ttg-breast, LR exhibits more stable ranking replicates than MLP and GNN
(diagonal of the heatmap). As expected, the more constant the classification
performance (Table 2), the more similar the ranking between replicates. Across
methods, the top-ranked genes may differ significantly. Classifier-based
methods, except GNN, tend to select similar genes, possibly reflecting linear
relationships with the class. LR with the embedded (weight) and post-hoc (IG)
scores generate distinct rankings, especially for LR+L2. As IG is also applied
to MLP and GNN, results across ML models are often more consistent when LR
(IG) is considered instead of LR (weight). Thus, in the following, experiments
are only conducted with IG. For the other datasets, qualitatively similar
results are observed, with even greater differences between the selected gene
sets obtained by ML models and statistical methods.

###  3.3 Explaining the predictions with a small set of genes

PGU100 - PGI$0$$20$$40$LR+L1LR+L2MLPGNNProportion of genes (%)(a) ttg-
breast$0$$20$$40$Proportion of genes (%)(b) BRCA-pam Figure 3: Impact of
progressive gene masking on the predictions of ML models (experiment 0). Genes
are masked by increasing (PGU) or decreasing order of importance (PGI) based
on the rankings $\boldsymbol{\phi}^{\text{IG}}$. For each data sample, PGU
calculates the percentage of well-ranked genes that should remain unmasked to
avoid disturbing a trained model. $100-$PGI estimates the percentage of well-
ranked genes that can be masked before disturbing the model. PGs are averaged
over all training samples correctly classified. Error bars are standard
deviations across replicates.

We investigate the minimum number of genes needed for a model to attain
performance comparable to using all genes. With LR+L1, the optimisation
process promotes sparsity, leading to many null parameters and the selection
of a small number of crucial genes, as shown in Table 3. The precise number
depends heavily on the choice of the hyperparameter $\lambda$. Still, the
overall strong performance of LR+L1 indicates the ability to classify
phenotypes effectively in these datasets with a limited number of genes.

Experiments are conducted to explore ML models. Experiment 0 is focused on the
ML-based rankings derived from IG. The number of genes used by a trained
classifier to make a decision is estimated by progressively masking the
lowest-ranked genes without re-training. Experiment 1 measures the performance
of models when they are trained from scratch using only the top-ranked genes
identified by IG and statistical methods. The goal is to assess whether the
top-ranked genes is sufficient for classification.

For experiment 0, in Fig. 3, the Prediction Gap on Unimportant features (PGU)
measures the minimum percentage of top-ranked genes that must remain unmasked
to avoid disturbing the model’s decision. For ttg-breast, PGUs for LR+L1 is
about 0.5% of genes (around 70 genes), which is less than the number of non-
zero parameters in this model (Table 3). PGUs for MLP, GNN and LR+L2 are 10%,
19% and 25% of genes respectively (around 1400, 2700 and 3600 genes). For
BRCA-pam, LR+L1 needs 0.6% of genes (around 80 genes), while other models have
to keep a larger number of genes, between 24% and 27% (3300 to 3800 genes).
Results on the other datasets are shown in Fig. S2.

For experiment 1, accuracies computed after re-training the models are in Fig.
4. For ttg-breast, the best 10 (resp. 100) genes are sufficient for MLP, LR+L1
and LR+L2 (resp. GNN) to saturate the classification performance. For BRCA-
pam, 500 genes are sufficient for all models. On both datasets, LRs and MLPs
obtain better performance than GNNs, which are disadvantaged because the graph
structure is highly perturbed. Note that better classification performance can
be obtained for MLP when trained on gene sets selected by other methods (Fig.
5). Results on the other datasets are in Figs. S3 and S4.

Thus, for a model trained on all genes, most genes but a very small sets of
the _best_ genes can be masked without significantly perturbing predictions.
Even smaller sets of the _best_ genes are sufficient to train models with high
performance.

###  3.4 Unimportant genes matter

$1$$10$$100$$1\,000$$60$$80$$100$(a) ttg-breastBalanced accuracyFeatures
keptImportantUnimportant$1$$10$$100$$1\,000$$20$$40$$60$$80$$100$Number of
selected features(b) BRCA-pamBalanced accuracyModelsLR+L1LR+L2MLPGNN Figure 4:
Classification performance shown for models trained on features identified as
important (full lines, experiment 1) or unimportant (dashed lines, experiment
2). Balanced accuracies are reported as a function of the number of kept
features for ttg-breast (a) and BRCA-pam (b) datasets using the specified
models. Error bars are std from 10 replicates.

This analysis investigates whether the genes selected in the previous section
exclusively contain relevant information. For that, we conduct experiment 0 by
masking the highest-ranked genes without re-training the models and experiment
2 for which models are trained using the $n_{L}$ lowest-ranked genes.

For experiment 0, in Fig. 3, the Prediction Gap on Important features (PGI)
measures the minimum percentage of worst-ranked genes that must remain
unmasked to avoid disturbing the model’s decision. Here, the orange dots
represent $100-$PGI values, indicating the highest percentage of well-ranked
genes that can be masked before perturbing the model. When the orange dots
surpass the blue dots (PGU), it suggests that all the genes selected as
important in the previous subsection can be masked without disturbing the
model. Thus, the remaining information is sufficiently redundant. For ttg-
breast, this occurs for MLP only. For LR+L2, masking a small proportion of the
identified important genes is enough to disrupt the model. For BRCA-pam, LR+L1
can mask all the identified important elements. In contrast, for MLP, LR+L2,
and GNN, masking a small proportion of the important genes is sufficient.
These results show that when a large number of genes is perturbed, the model
is disrupted. However, when the model is based on a small number of genes,
these genes are not necessarily the only ones containing relevant information.
Results on the other datasets are in Fig. S2.

For experiment 2, with re-training, results are in Fig. 4 (dashed lines).
Using the 100 lowest-ranked genes is sufficient to achieve optimal balanced
accuracy on ttg-breast. For BRCA-pam, the gaps between the full lines (best-
ranked genes) and the dashed lines (worst-ranked genes) are small. Results on
the other datasets are in Fig. S3.

Similar classification performance can be achieved by keeping a relatively
small set of lower-ranked genes, showing there is no unique set of informative
genes.

$1$$10$$100$$1\,000$$60$$80$$100$(a) ttg-breastBalanced
accuracyFiltersDESeq2MIVARPCAEdgeR$1$$10$$100$$1\,000$$20$$40$$60$$80$$100$Number
of selected features(b) BRCA-pamBalanced accuracyIGLR+L1LR+L2MLPGNN Figure 5:
Classification performance of a MLP trained on sets of features identified as
important by various methods as indicated. The representation is coded as in
Fig. 4.

((a)) ttg-breast

((b)) BRCA-pam

Figure 6: Top 10 over-represented gene sets identified by the different
methods (blue dots). When one of these gene sets appear among the top 100
over-represented sets of another method, it is flagged with a light green dot.

###  3.5 Stability across gene sets

The top-ranked genes identified by the different methods are compared with the
content of gene sets known for their biological properties using an over-
representation analysis. For each method, the top 10 over-represented gene
sets are displayed in Fig. 6 (blue dots). When these gene sets appear among
the top 100 significant gene sets identified by another method, they are
flagged with a light green dot. Interestingly, different sets associated with
(breast) cancer emerge across different methods. However, the overlap between
the methods can be small. See Fig. S5 for the other datasets.

##  4 Discussion

Many ML models have emerged to classify phenotypes from gene expression data,
with the goal of profiling genes discriminative for the phenotypes. Model
explanations are often presented as a list of genes ranked by importance. Our
objective is to evaluate the relevance of explanations generated by ML models,
examining whether they contribute novel perspectives to elucidate complex
biological pathways.

To address this question, experiments are carried out on human tissue samples
coming from TCGA, TARGET and GTEx databases. Three ML models, LRs, MLPs and
GNNs, are trained to discern a pathological state from other pathological or
healthy states. Gene rankings are determined using the integrated gradient
method (IG), a neural network-tailored explainability technique, and are
compared with 3 commonly-used statistical methods: DESeq2, EdgeR and mutual
information (MI). Additionally, an over-representation analysis is conducted
to explore the biological relevance of the top-ranked genes.

Several key insights emerge. First, the classifiers consistently achieve a
balanced accuracy exceeding 95% in the majority of cases, with best
performance being obtained by LR. Additionally, as cross-validation
experiments favoured shallow MLP comprising 1 or 2 layers, this suggests
limitations in leveraging complex relationships within the data. Second, the
classifiers can maintain the same performance using only on a small set of
top-ranked genes for training. These genes are potential biomarkers for the
studied pathologies. Yet, substantial differences appear in the genes
identified across the different methods, highlighting the influence of
inherent bias in each gene selection method.

However, good classification performance is also obtained with limited sets of
ranked lower, revealing both dispersion and redundancy of information in the
gene space. This dispersion is easily measurable for the binary classification
tasks with a Welch’s unequal variance t-test corrected for multiple testing.
This statistical test quantifies the difference between the means of the
distributions for each class, accounting for their average variances. Taking
the ttg-breast dataset as an example, around 90% of the genes have an adjusted
p-value below 5%. Most genes carry enough information to distinguish the
classes. Hence, in trained models, up to 20 % of the most important genes
($100-PGI\simeq 20\%$) can be masked without affecting classification
performance. Similarly, re-training models on lower-ranked genes achieves good
performance.

To gain insights into the genes preferentially used by ML methods, looking at
the correlation of rankings with the t-statistic is informative. On the ttg-
breast dataset, scores from DESeq2, EdgeR and MI are highly correlated with
the t-statistic, with Spearman correlations around $0.9$. LR and MLP exhibit
lower correlation around $0.45$. However, on a more complex dataset such as
ttg-all, characterised by more heterogeneous data samples, correlations of the
ML-based rankings with the t-statistic decrease. Surprisingly, GNN
consistently exhibits low correlations, notably for ttg-all where the
t-statistic distribution of the top 100 genes selected by IG is similar to
those of genes chosen randomly. This highly questions the interpretability of
genes by the GNN (IG) method.

In terms of classification performance, statistical methods looking for
differentially expressed genes prove to be even better than ML methods in
selecting a minimal set of informative genes (Fig. 5). A hypothesis is that
during training, the MLP performance rapidly converges to 100% of accuracy,
making the use of the most most differentially expressed genes unnecessary.

From a biological point of view, the concept of explainability, as defined
within the ML community, does not immediately provide a novel perspective to
unravel complex biological pathways. Given the dispersed and redundant nature
of relevant information, transitioning from individual genes to an exploration
of cellular processes holds promise for enhancing our understanding of
biological phenomena. Here, the over-representation analysis of the top 100
genes identified by the different methods highlights certain cellular
processes, occasionally related to the studied pathologies. While some of them
may be shared among methods, there is a substantial variability. The different
genes identified by the methods may prove complementary, underscoring distinct
processes leading to a phenotype. Further investigation of functional sets of
genes, using for instance mechanistic experiments, could help to confirm this
hypothesis.

In summary, this study provides valuable insights for researchers seeking
biologically relevant molecular signatures of pathologies using ML. The
investigation of functional gene sets is a promising research direction.

##  5 Competing interests

No competing interest is declared.

##  6 Author contributions statement

MyB, JA, BA, PB designed the work. MyB, AH, MaB mined the bioinformatics
literature. MyB retrieved the datasets and conceived the experiments. MyB
prepared the first draft of the manuscript, the figures and output tables. All
authors analysed the result, revised, read and approved the final manuscript.

##  7 Acknowledgments

Grants: CHIST-ERA-19-XAI-006 and GRAPHNEX ANR-21-CHR4-0009. We thank the
Centre Blaise Pascal’s IT test platform (ENS de Lyon, France) for the ML
facilities quemener2013sidus .

## References

  * [1] MD Robinson, DJ McCarthy, et al.  edgeR: a bioconductor package for differential expression analysis of digital gene expression data.  Bioinformatics, 26(1):139–140, 2010.
  * [2] MI Love, W Huber, et al.  Moderated estimation of fold change and dispersion for rna-seq data with DESeq2.  Genome biology, 15(12):1–21, 2014.
  * [3] I Goodfellow, Y Bengio, and Aaron Courville.  Deep Learning.  MIT Press, 2016.  http://www.deeplearningbook.org.
  * [4] TR Golub, DK Slonim, et al.  Molecular classification of cancer: class discovery and class prediction by gene expression monitoring.  Science, 286(5439):531–537, 1999.
  * [5] KA Hoadley, C Yau, et al.  Cell-of-origin patterns dominate the molecular classification of 10,000 tumors from 33 types of cancer.  Cell, 173(2):291–304, 2018.
  * [6] JS Parker, M Mullins, et al.  Supervised risk predictor of breast cancer based on intrinsic subtypes.  Journal of clinical oncology, 27(8):1160, 2009.
  * [7] D Leng, L Zheng, et al.  A benchmark study of deep learning-based multi-omics data fusion methods for cancer.  Genome biology, 23(1):1–32, 2022.
  * [8] B Hanczar, V Bourgeais, et al.  Assessment of deep learning and transfer learning for cancer prediction based on gene expression data.  BMC bioinformatics, 23(1):262, 2022.
  * [9] Y Li, K Kang, et al.  A comprehensive genomic pan-cancer classification using the cancer genome atlas gene expression data.  BMC genomics, 18:1–13, 2017.
  * [10] E Jacquet, F Chuffart, et al.  Aberrant activation of five embryonic stem cell-specific genes robustly predicts a high risk of relapse in breast cancers.  BMC genomics, 24(1):463, 2023.
  * [11] TNR Rohimat, F Nhita, et al.  Implementation of genetic algorithm-support vector machine on gene expression data in identification of non-small cell lung cancer in nonsmoking female.  In IC2IE, pages 361–366, 2022.
  * [12] T Ahn, T Goo, et al.  Deep learning-based identification of cancer or normal tissue using gene expression data.  In BIBM, pages 1748–1752. IEEE, 2018.
  * [13] R Ramirez, Y Chiu, et al.  Classification of cancer types using graph convolutional neural networks.  Frontiers in physics, 8:203, 2020.
  * [14] N Mahendran, PM Durai Raj Vincent, et al.  Machine learning based computational gene selection models: a survey, performance evaluation, open issues, and future research directions.  Frontiers in genetics, 11:603808, 2020.
  * [15] F Alharbi and A Vakanski.  Machine learning methods for cancer classification using gene expression data: A review.  Bioengineering, 10(2):173, 2023.
  * [16] C Molnar.  Interpretable Machine Learning.  2022\.
  * [17] P Khatri, M Sirota, et al.  Ten years of pathway analysis: current approaches and outstanding challenges.  PLoS computational biology, 8(2):e1002375, 2012.
  * [18] M Sundararajan, A Taly, et al.  Axiomatic attribution for deep networks.  In ICML, pages 3319–3328. PMLR, 2017.
  * [19] J Lonsdale, J Thomas, et al.  The genotype-tissue expression (gtex) project.  Nature genetics, 45(6):580–585, 2013.
  * [20] CM Bishop.  Pattern recognition and machine learning.  Springer, 2006.
  * [21] J Zhou, G Cui, et al.  Graph neural networks: A review of methods and applications.  AI open, 1:57–81, 2020.
  * [22] TN Kipf and M Welling.  Semi-supervised classification with graph convolutional networks.  In ICLR, 2017.
  * [23] IS Dhillon, Y Guan, et al.  Weighted graph cuts without eigenvectors a multilevel approach.  IEEE TPAMI, 29(11):1944–1957, 2007.
  * [24] IT Jolliffe.  Principal Component Analysis.  Springer, New York, NY, second edition, 2002.
  * [25] BC Ross.  Mutual information between discrete and continuous data sets.  PloS one, 9(2):e87357, 2014.
  * [26] J Costa-Silva, D Domingues, et al.  Rna-seq differential expression analysis: An extended review and a software tool.  PloS one, 12(12):e0190152, 2017.
  * [27] SM Lundberg and S Lee.  A unified approach to interpreting model predictions.  NeurIPS, 30, 2017.
  * [28] C Agarwal, S Krishna, et al.  Openxai: Towards a transparent evaluation of model explanations.  NeurIPS, 35:15784–15799, 2022.
  * [29] A Liberzon, A Subramanian, et al.  Molecular signatures database (MSigDB) 3.0.  Bioinformatics, 27(12):1739–1740, 2011.
  * [30] A Liberzon, C Birger, et al.  The molecular signatures database hallmark gene set collection.  Cell systems, 1(6):417–425, 2015.
  * [31] RL Grossman, AP Heath, et al.  Toward a shared vision for cancer genomic data.  New England Journal of Medicine, 375(12):1109–1112, 2016.
  * [32] MJ Goldman, B Craft, et al.  Visualizing and interpreting cancer genomics data via the xena platform.  Nature biotechnology, 38(6):675–678, 2020.
  * [33] K Tomczak, P Czerwińska, et al.  Review the cancer genome atlas (TCGA): an immeasurable source of knowledge.  Contemporary Oncology, 2015(1):68–77, 2015.
  * [34] GF Gao, JS Parker, et al.  Before and after: comparison of legacy and harmonized tcga genomic data commons’ data.  Cell systems, 9(1):24–34, 2019.
  * [35] J Vivian, AA Rao, et al.  Toil enables reproducible, open source, big biomedical data analyses.  Nature biotechnology, 35(4):314–316, 2017.
  * [36] S Rhee, S Seo, et al.  Hybrid approach of relation network and localized graph convolutional filtering for breast cancer subtype classification.  In IJCAI, pages 3527–3534, 2018.
  * [37] JM Choi and H Chae.  moBRCA-net: a breast cancer subtype classification framework based on multi-omics attention neural networks.  BMC bioinformatics, 24(1):1–15, 2023.
  * [38] F Pedregosa, G Varoquaux, et al.  Scikit-learn: Machine learning in Python.  JMLR, 12:2825–2830, 2011.
  * [39] A Paszke, S Gross, et al.  Pytorch: An imperative style, high-performance deep learning library.  NeurIPS, 32, 2019.
  * [40] DP Kingma and J Ba.  Adam: A method for stochastic optimization.  In Yoshua Bengio and Yann LeCun, editors, ICLR, 2015.
  * [41] DJB Clarke, M Jeon, et al.  Appyters: Turning jupyter notebooks into data-driven web apps.  Patterns, 2(3), 2021.
  * [42] N Kokhlikyan, V Miglani, et al.  Captum: A unified and generic model interpretability library for pytorch.  arXiv preprint arXiv:2009.07896, 2020.
  * [43] E Quemener and M Corvellec.  SIDUS—the solution for extreme deduplication of an operating system.  Linux Journal, 2013(235):3, 2013.

##  8 Supplementary material

Datasets | LR | MLP | GNN
---|---|---|---
$\lambda$ (L1) |  $\lambda$ (L2) | Lay. | Feat. | Lay. | Feat. | $k$
PanCan | 1 | 0.1 | 1 | 20 | 1 | 2 | 2
BRCA | 0.1 | 0.01 | 1 | 20 | 1 | 2 | 2
BRCA-pam | 0.1 | 1 | 1 | 20 | 1 | 2 | 10
ttg-breast | 0.1 | 0.1 | 1 | 20 | 1 | 1 | 10
ttg-all | 1 | 0.1 | 2 | 40 | 1 | 2 | 2

Table S1: Hyperparameters selected by grid search: the weight $\lambda$ for
the LR regularisation term (for trade-off with the cross-entropy loss)
($\lambda$ chosen as 0.01, 0.1, 1 or 10); the number of MLP layers (1 or 2),
their number of hidden features (10, 20, 40 or 80); the number of GNN layers
(1, 3 or 5), their number of hidden features (1 or 2), and $k$ setting the
number of graph edges to $k\times G$ ($k$ is 2 or 10).

Dataset | LR+L1 | LR+L2 | MLP | GNN
---|---|---|---|---
PanCan | $96.6$ | $96.2$ | $96.0\pm 0.2$ | $94.5\pm 0.2$
BRCA | $99.4$ | $97.3$ | $99.3\pm 0.1$ | $99.2\pm 0.1$
BRCA-pam | $91.8$ | $89.5\pm 0.2$ | $87.7\pm 1.3$ | $87.5\pm 1.3$
ttg-breast | $99.5$ | $98.7\pm 0.1$ | $99.3\pm 0.2$ | $99.2\pm 0.1$
ttg-all | $99.5$ | $99.5$ | $99.6$ | $99.4\pm 0.1$

Table S2: Classification performance measured by accuracy (%). Standard
deviations are computed from 10 replicates; not reported when below 0.05.

Dataset | LR+L1 | LR+L2 | MLP | GNN
---|---|---|---|---
PanCan | 58253 | 42564 | 74 | 404
BRCA | 222 | 175 | 6 | 8
BRCA-pam | 673 | 496 | 7 | 11
ttg-breast | 262 | 207 | 10 | 12
ttg-all | 6127 | 2703 | 136 | 150

Table S3: Average training duration (s) for each model.

((a)) BRCA

((b)) PanCan

((c)) ttg-all

Figure S1: Heatmaps showing the percentage of common genes among the top 10
(lower) and top 100 (upper + diagonal) genes selected by each method.
$0$$20$$40$LR+L1LR+L2MLPGNNProportion of genes (%)(a)
BRCA$0$$20$$40$$60$$80$Proportion of genes (%)(b) ttg-
all$0$$50$$100$Proportion of genes (%)(c) PanCanPGU100 - PGI Figure S2: Impact
of progressive gene masking on the predictions of ML models. Genes are masked
by increasing (PGU) or decreasing importance (PGI) using
$\boldsymbol{\phi}^{\text{IG}}$ scores. PGs are averaged over all training
samples correctly classified without masking, with error bars representing
standard deviations.  $1$$10$$100$$1\,000$$60$$80$$100$Number of selected
features(a) BRCABalanced accuracy$1$$10$$100$$1\,000$$60$$80$$100$Number of
selected features(b) ttg-all$1$$10$$100$$1\,000$$0$$20$$40$$60$$80$$100$Number
of selected features(c) PanCanBalanced accuracyFeatures
keptModelsImportantLR+L1UnimportantLR+L2MLPGNN Figure S3: Classification
performance shown for models trained on features identified as important or
unimportant for each model. Balanced accuracies are reported as a function of
the number of features kept using the specified models. Error bars are std
from 10 replicates.  $1$$10$$100$$1\,000$$60$$80$$100$Number of selected
features(a) BRCABalanced accuracy$1$$10$$100$$1\,000$$60$$80$$100$Number of
selected features(b) ttg-all$1$$10$$100$$1\,000$$0$$20$$40$$60$$80$$100$Number
of selected features(c) PanCanBalanced accuracyFiltersPost-
hocDESeq2LR+L1MILR+L2VARMLPPCAGNNEdgeR Figure S4: Classification performance
of a MLP trained on sets of features identified as important by various
methods.

((a)) BRCA

((b)) ttg-all

((c)) PanCan

Figure S5: Top 10 over-represented gene sets identified by the different
methods (blue dots). When one of these gene sets appear among the top 100
over-represented sets of another method, it is flagged with a light green dot.