Conformal novelty detection for multiple metabolic networks

Abstract

Background

Graphical representations are useful to model complex data in general and biological interactions in particular. Our main motivation is the comparison of metabolic networks in the wider context of developing noninvasive accurate diagnostic tools. However, comparison and classification of graphs is still extremely challenging, although a number of highly efficient methods such as graph neural networks were developed in the recent decade. Important aspects are still lacking in graph classification: interpretability and guarantees on classification quality, i.e., control of the risk level or false discovery rate control.

Results

In our contribution, we introduce a statistically sound approach to control the false discovery rate in a classification task for graphs in a semi-supervised setting. Our procedure identifies novelties in a dataset, where a graph is considered to be a novelty when its topology is significantly different from those in the reference class. It is noteworthy that the procedure is a conformal prediction approach, which does not make any distributional assumptions on the data and that can be seen as a wrapper around traditional machine learning models, so that it takes full advantage of existing methods. The performance of the proposed method is assessed on several standard benchmarks. It is also adapted and applied to the difficult task of classifying metabolic networks, where each graph is a representation of all metabolic reactions of a bacterium and to real task from a cancer data repository.

Conclusions

Our approach efficiently controls — in highly complex data — the false discovery rate, while maximizing the true discovery rate to get the most reasonable predictive performance. This contribution is focused on confident classification of complex data, what can be further used to explore complex human pathologies and their mechanisms.

Background

With the rise of new sequencing and high-throughput technologies, new data in form of metabolic networks are more and more available to support the study of human pathologies. These datasets are huge and of complex structure requiring the application of appropriate machine learning models. In particular, the interest to apply predictive models to metabolic information is extremely high in metabolic diseases tasks, such as prediction of obesity and diabetes. These pathologies result from a certain disability of a cell to breakdown or produce some essential substrates. As a result, if an enzyme in one reaction is broken, it may influence subsequent reactions, leading to enormous cascading damages [1, 2].

A metabolic network represents a complete set of metabolic and physical processes describing physiological and biochemical properties of a living cell [3]. Moreover, modern large metabolic databases such as KEGG [4] make it possible to access genomic, enzymatic and metabolic information and to reconstruct interactions. The representation of a graph is an important question. Adjacency matrices are a natural representation of metabolic networks, and the Graph Neural Networks (GNN) which we consider below, cope very reasonably with the complex data analysis. Extraction of the topological properties of graphs is another efficient solution, since it reduces the data dimension and often leads to the state-of-the-art performance [5]. So, strong correlations between phenotypical traits of organisms and the topology of metabolic networks were reported [6, 7], underlining the importance to study metabolic networks and their topology. An efficient classification of patient-specific metabolic networks from a breast cancer study using topological properties was reported by [5], illustrating a strong motivation to make use of metabolic networks in medical tasks.

From a mathematical viewpoint, a metabolic network can be represented by a graph composed of nodes and edges, which connect the nodes. The metabolites and enzymes are the nodes of the graph. Each reaction substrate is linked to the catalysing enzyme and each enzyme is connected to products of the chemical reaction  [8]. Modern statistical and machine learning methods can be applied to such data to get more insights in the functioning of living organisms [9].

While the overwhelming majority of existing results concern the analysis of a single network, here we are particularly interested in the classification and comparison of multiple metabolic networks, since it is a step forward to the development of non invasive accurate diagnostic tools. Our specific goal is novelty detection, sometimes also called outlier detection. This amounts to compare metabolic networks to a reference and to decide which of the networks are significantly different from the reference. Detecting metabolic networks whose structure or topology is inherently different from the structure of a set of default or nominal graphs is crucial for the identification of anormal cells and for gaining new insights into the functionalities of different metabolisms.

Network comparison is an inherently involved problem due to the complex structure of graphs. Generally, dimensionality reduction methods are used which provide a graph embedding for every network. This can be achieved by traditional principal component analysis [10, 11] or more recently with graph neural networks [12,13,14,15]. In medical and pre-clinical research, novelty detection based on graph embeddings is generally done in a manual way, which is both time consuming for human experts and highly subjective, since a human might take a biased decision instead of using an objective criterion. Moreover, such novelty detection comes without any guarantee on the quality of the results.

In many applications, it is vital not only to make accurate predictions but also to quantify the accuracy and provide explanations on the learned model. Here, we follow the definition proposed by [16, 17] that a novelty is any newly acquired structure or property that permits the assumption of a new function. With the hope that studies about metabolic alterations in organisms help to decipher mechanisms of pathologies, [18] discusses novelty detection and tests several classification methods to detect novel metabolites of a breast cancer-tissue sample that can appear during dynamic evolution of biological cells. It was also recently reported that metabolic changes, i.e., novelties, can help in diagnostic and therapeutic investigations; a repository of tumor metabolic networks derived from Genome-Scale Metabolic Models, called TumorMet [19], is publicly available.

More precisely, any novelty detection method may make mistakes either by declaring observations as novelties while they are not, or by not recognizing a new observation as a novelty. Depending on the difficulty of the underlying problem, that is, whether novelties are completely different from the reference or still share some similarities, the number of errors for a given method may vary greatly. Traditional machine learning methods do not provide any information on the quality of its results. However, recently, new procedures have been developed that come with a statistical guarantee that the set of detected novelties contains at most a given percentage of falsely detected novelties, while keeping the number of identified novelties as large as possible. One of such approaches is conformal prediction [20, 21] first introduced in classification and regression settings. Figure 1 illustrates the novelty detection task and the possible error types that a procedure can make.

A major advantage of conformal prediction is that it does not make any assumptions on the type of distribution of the observations, which is important for applications where distribution assumptions on the data are hard to verify and/or rarely satisfied. Moreover, the reference or nominal distribution is not assumed to be known, but it is sufficient to dispose of a sample from the reference distribution, that is, a semi-supervised setting is considered. Here, we use the notion of the semi-supervised scenario introduced by [22]. Such a general framework makes conformal prediction highly attractive for uncertainty quantification on complex structures such as networks. Moreover, conformal prediction is known to provide non-asymptotic and distribution-free coverage guarantees for various tasks [23, 24]. Finally, it is extremely noteworthy that conformal prediction is a wrapper around traditional machine learning models, that is, it can be directly applied to the output of existing methods. As such, conformal prediction takes full advantage of existing high-performance machine learning algorithms.

Fig. 1

Illustration of the data sets, the novelty detection task and possible errors. All networks in the reference set are composed of two communities, while the test set also contains networks with one community. The novelty detection method identifies most of the novelties (in red) correctly, but also falsely declares a reference network as a novelty

In a recent series of works, conformal procedures for the novelty detection task have been developed [22, 25,26,27,28,29]. The general idea is to, first, learn a non-conformity score for all new observations using some existing machine learning algorithm. Then a comparison of these scores to the scores of the reference observations provides the final set of detected novelties. For the selection of the final set of novelties, tools from multiple testing are used, where finite-sample guarantees on the error rate can be obtained. The procedure AdaDetect proposed by [30] is the most powerful approach up to date, which is based on a particularly efficient way of learning the non-conformity scores, and it is appropriate for a huge variety of settings.

Novelty, or outlier, detection in graphs is particularly challenging due to the complex structure of the data. Some very recent attempts to extend conformal prediction to graphs concern either link prediction [31] or the prediction of node labels [32, 33] in a single graph. A recently developed procedure for the problem of novelty detection in a collection of networks is given in [34], where a simple outlier detection method is proposed; this approach is applied to graphs in neuroscience and relies on a hierarchical generalized linear model, but not on conformal inference.

Our aim here is to extend AdaDetect to the specific task of novelty detection in a collection of metabolic networks. Our goal is to show its usefulness in practice and illustrate the gain of new insights on the cell metabolism. Our contribution is multi-fold:

• We propose a statistically sound method that identifies the networks in a data set which are significantly different from a set of provided reference observations and that controls the corresponding false discovery rate. Our approach applies to any type of complex networks, not necessarily to metabolic networks. To our best knowledge, we are the first ones to propose conformal outlier detection procedure for a collection of graphs.

• We discuss theoretical guarantees of our method, and the ability of the proposed procedure to outperform some existing state-of-the-art baseline methods.

• We show the generalizing performance of the newly introduced model on several benchmarks.

• We validate our method on a real data set of bacteria, where each bacterium is represented by its metabolic network. We also apply the proposed approach to two tasks – Lung and Kidney cancer – of the TumorMet repository.

Note that both novelty and outlier detection are methods to identify anomalies. However, traditionally, in an outlier detection task a dataset contains outliers which should be detected. Therefore, the outlier detection is an unsupervised anomaly detection approach. In a novelty detection task, a learning algorithm has access to the nominal class of observations, and the novelties are not observed. Such a novelty detection scenario is called semi-supervised, since only one class of labeled data is available. More details on the semi-supervised novelty detection framework – and our contribution is in this context – can be found in [35].

Method: AdaDetect for Graphs

Setting and notations

We consider the general setting with observed networks denoted by $G=(A,X)$, where A is the adjacency matrix of the network and X is a matrix of node covariates (if available). Networks may be directed or undirected, binary or valued, share the same nodes over all networks or not, may have varying number of nodes from one network to the other or include node covariates.

In the semi-supervised framework, two sets of networks are observed. First, ${{\,\mathrm{\mathcal {G}_{\text {ref}}}\,}}=\{G_i, i\in \llbracket 1,n\rrbracket \}$ is a set of networks having the standard or normal behavior, referred to as the reference sample. Here $\llbracket a,b\rrbracket$ denotes the set of integers from a to b. These networks are assumed to be i.i.d. realizations of some distribution $P_\text {ref}$, which is unknown to the user. That is, $G_i \sim P_\text {ref}, i\in \llbracket 1,n\rrbracket$. The second set of observed networks denoted by ${{\,\mathrm{\mathcal {G}_{\text {test}}}\,}}=\{G_{i}, i\in \llbracket n+1,n+m\rrbracket \}$ is the test sample, where the observed networks are assumed to be independent, not observed during training. The task is to decide which of them are novelties, that is, which networks do not come from the reference distribution $P_\text {ref}$. To put it differently, the aim is to discover the set $\mathcal {I}_{\textrm{nov}} = \{ i\in \llbracket n+1, n+m\rrbracket , G_i \not \sim P_\text {ref} \}$, which is the set of indices of the novelties. Furthermore, denote ${{\,\mathrm{\mathcal {I}_{\text {test}}}\,}}=\llbracket n+1,n+m\rrbracket$ the set of indices of networks in ${{\,\mathrm{\mathcal {G}_{\text {test}}}\,}}$ and $\mathcal {I}_{\textrm{ref}}={{\,\mathrm{\mathcal {I}_{\text {test}}}\,}}\backslash \mathcal {I}_{\textrm{nov}}$ the set of indices i of networks in ${{\,\mathrm{\mathcal {G}_{\text {test}}}\,}}$ from the reference distribution, that is $G_i\sim P_\text {ref}$.

Now a novelty detection procedure is a (measurable) function R that returns a subset of ${{\,\mathrm{\mathcal {I}_{\text {test}}}\,}}$ corresponding to the indices of the networks declared as novelties. The false discovery rate (FDR), that is the proportion of falsely declared novelties, and the true discovery rate (TDR), that is the proportion of correctly identified novelties, are defined as

$$FDR(R) = {\mathbb{E}} \left[ {\frac{{|{\mathcal{I}}_{{ref}} \cap R|}}{{1 \vee |R|}}} \right], TDR(R) = {\mathbb{E}} \left[ {\frac{{|{\mathcal{I}}_{{nov}} \cap R|}}{{1 \vee |{\mathcal{I}}_{{nov}} |}}} \right]$$

Our goal is to find a procedure that controls the FDR at a prescribed level α, that is ${{\text{FDR}}}(R) \le \alpha$, and whose power or TDR is as large as possible.

AdaDetect

The general idea of conformal novelty detection [20, 22, 25, 36] is to compare the non-conformity score of a test observation to the scores of the reference observations to decide whether the observation is a novelty or not. In the seminal work of [25], the FDR control is shown to be guaranteed when proceeding as follows: the reference sample ${{\,\mathrm{\mathcal {G}_{\text {ref}}}\,}}$ is split into two parts and the score is trained using one of the parts. The power of the procedure, that is the number of correctly detected novelties, depends on the quality of the scores. In this line, [27] recently achieved an important improvement by training the score not only on the reference sample ${{\,\mathrm{\mathcal {G}_{\text {ref}}}}}$ with a one-class classification algorithm, but training the score using a two-class classification method including the test sample ${{\,\mathrm{\mathcal {G}_{\text {test}}}\,}}$. This procedure is called AdaDetect and yields a significant gain in power, while still controlling the FDR. In this work, we work out how to successfully apply AdaDetect to the task of novelty detection in network data.

Algorithm 1

AdaDetect for networks

Fig. 2

Choice of threshold $\delta$ for the novelty detection procedure $R_\delta$. Illustration of proportion $p_\delta$ of falsely declared novelties for three different thresholds

In detail, in AdaDetect the non-conformity score is a classifier trained on the following problem. First, the reference sample ${{\,\mathrm{\mathcal {G}_{\text {ref}}}\,}}$ is split into two subsets, say ${{\,\mathrm{\mathcal {G}_{\text {train}}}\,}}$ and ${{\,\mathrm{\mathcal {G}_{\text {cal}}}\,}}$. The networks in ${{\,\mathrm{\mathcal {G}_{\text {train}}}\,}}$ are labeled as “0” and the elements of ${{\,\mathrm{\mathcal {G}_{\text {cal}}}\,}}\cup {{\,\mathrm{\mathcal {G}_{\text {test}}}\,}}$ as “1”. That is, ${{\,\mathrm{\mathcal {G}_{\text {train}}}\,}}$ is a pure class, where all observations come from the reference distribution $P_\text {ref}$, while ${{\,\mathrm{\mathcal {G}_{\text {cal}}}\,}}\cup {{\,\mathrm{\mathcal {G}_{\text {test}}}\,}}$ is mixed, containing observations from the reference distribution $P_\text {ref}$ as well as novelties. A classifier can be trained using any machine learning classifier for graphs. The classifier function, say g, gives the non-conformity score, which generally corresponds to the probability of a network to belong to class “1”.

Now denote by ${{\,\mathrm{\mathcal {I}_{\text {cal}}}\,}}\subset {{\,\mathrm{\mathcal {I}_{\textrm{ref}}}\,}}$ the set of indices for which $G_i \in {{\,\mathrm{\mathcal {G}_{\text {cal}}}\,}}$. For every $i\in {{\,\mathrm{\mathcal {I}_{\text {cal}}}\,}}\cup {{\,\mathrm{\mathcal {I}_{\text {test}}}\,}}$, compute the non-conformity score $S_i = g(G_i)$. Now it is reasonable to declare all networks that have a large score, or more precisely, whose score is larger than some threshold $\delta$, as novelties. That is, we consider the novelty selection procedure $R_\delta$ defined as

$$\begin{aligned} R_\delta =\{j\in {{\,\mathrm{\mathcal {I}_{\text {test}}}\,}}, S_j\ge \delta \}. \end{aligned}$$

As the choice of the threshold is crucial for the control of the FDR, we introduce the quantity $p_\delta$ defined as the proportion of reference observations declared as novelties by procedure $R_\delta$:

$$\begin{aligned} p_\delta = \frac{ |\{ i\in {{\,\mathrm{\mathcal {I}_{\text {cal}}}\,}}: S_i > \delta \}| +1 }{|{{\,\mathrm{\mathcal {I}_{\text {cal}}}\,}}|+1}. \end{aligned}$$

That is, $p_\delta$ is an approximation of the FDR of procedure $R_\delta$. Thus, for our procedure we choose the smallest threshold $\delta$ such that the corresponding proportion $p_\delta$ is lower than $\alpha$, which is the desired risk level. The justification of this approach is that it is impossible for any classifier to learn to distinguish reference observations in class “0” from reference observations in class “1”. Thus, the distribution of the scores $\{S_i,i\in {{\,\mathrm{\mathcal {I}_{\text {cal}}}\,}}\}$ equals the distribution of the score under $P_\text {ref}$. Hence, $\{S_i,i\in {{\,\mathrm{\mathcal {I}_{\text {cal}}}\,}}\}$ is appropriate for a comparison of the scores of the test observations for fixing the threshold $\delta$ that yields the FDR control.

The procedure is summarized in Algorithm 1. Figure 2 schematically shows the proportion of falsely declared novelties for different thresholds. [27] illustrate that AdaDetect is more powerful than state-of-the art procedures, that is, it detects more novelties than other methods.

Note that while the FDR control holds regardless of sample sizes, in practice, the sample size of the calibration set ${{\,\mathrm{\mathcal {G}_{\text {cal}}}\,}}$ must be large enough to ensure a good power [22, 25, 27]. However, increasing $\vert {{\,\mathrm{\mathcal {G}_{\text {cal}}}\,}}\vert$ and hence the proportion of references networks in the mixed set ${{\,\mathrm{\mathcal {G}_{\text {cal}}}\,}}\cup {{\,\mathrm{\mathcal {G}_{\text {test}}}\,}}$ degrades the quality of the scores learned by the classifier in AdaDetect. Based on power results from [22, 27], we recommend to choose $\vert {{\,\mathrm{\mathcal {G}_{\text {cal}}}\,}}\vert$ to be of the same order as the test sample size $|{{\,\mathrm{\mathcal {G}_{\text {ref}}}\,}}|$. Ideally, ${\mathcal {G}}_{cal}$ has to be as large as possible. However, we defined ${\mathcal {G}}_{cal} = {\mathcal {G}}_{train} = n/2$, where n is the limited number of available observations. If we increase ${\mathcal {G}}_{cal}$, our ${\mathcal {G}}_{train}$ will be smaller what will lead to a worse performance.

Machine learning algorithms for graph classification

Graph classification has garnered significant interest in recent years with many new methods, especially in the field of deep learning [37]. A particular feature of graphs is that in general no natural order of the vertices exits, but that classifiers are required to be invariant to the reordering of the nodes. This is one of the main obstacles to extending classical ML methods to graph data. Hence, as networks are complex data objects, their comparison is a challenging task and different methods handle this question in widely different ways. The choice of an appropriate approach also depends on the characteristics of the data at hand such as the availability of node features, the absence or presence of node correspondence, or whether the networks are directed or not. In this section, we present the general approaches to graph classification and we discuss their principal properties. The two main approaches are graph kernels [38] and graph neural networks (GNNs) [37, 39,40,41,42].

Graph kernels are the historically dominant approach for graph classification. A graph kernel is a deterministic function defining a similarity measure between a pair of networks, that can be combined with a support vector machine (SVM) classifier for supervised learning. The most popular graph kernel is the Weisfeiler-Lehman (WL) kernel [43], which is suited for networks with discrete node attributes. It is based on the 1-WL or color refinement algorithm, which proceeds as follows: First, for every network a graph embedding is computed according to some message passing mechanism. In detail, for every node its label (or color) is aggregated with the labels of its neighbors yielding a fingerprint. Then all nodes with identical fingerprint are assigned a new common node label (or color).

When this procedure is iterated, say K times, this results in a node embedding describing the structure of the K-hop neighborhood of every node, where nodes with identical K-hop neighborhoods share the same label. In other words, the node embeddings define a node clustering. For the WL subtree kernel [43] the clusterings obtained at all iterations are used to build a graph embedding with multiple resolutions. Finally, the WL kernel of two graphs is defined as the inner product of their graph embeddings. The WL algorithm gives rise to the most powerful existing test to decide whether two graphs are isomorphic, that is, whether one of the graphs can be obtained from the other by a permutation of the nodes [44].

Many other graph kernels have been proposed in the literature, see [38] for a complete review, some of them take into account discrete or continuous node attributes. A general main drawback of graph kernels is the computational burden that comes with computing the kernel function for every pair of networks in the training sample. The running time is ${\mathcal {O}}(NKn_{\max }+N^2Km_{\max })$, where $n_{\max }$ and $m_{\max }$ are the maximum number of vertices and edges in a collection of N graphs. In our case $N=\vert {{\,\mathrm{\mathcal {G}_{\text {cal}}}\,}}\vert + \vert {{\,\mathrm{\mathcal {G}_{\text {test}}}\,}}\vert$.

Graph neural networks (GNNs) are recent approaches for graph-based learning, that aim to scale to larger datasets than graph kernels by generalizing neural networks (NNs) to graph-structured data. Specifically, GNNs are permutation-invariant functions that produce a vector representation of each node in a network, using a combination of linear and non-linear operations. This vector representation is based on a neighborhood aggregation scheme, where, given an initial representation of the nodes, node representations are updated by taking into account their neighbors (e.g. computing the sum, mean or max of the representations). In that view, GNNs can be seen as a neural network version of the 1-WL algorithm. As GNNs produce a matrix of node embeddings, they can be combined with NN layers for node classification [45] or link prediction [46]. Moreover, graph classification can be performed by collapsing the node representations of a network into a single vector representation, and by feeding this graph representation into NN layers for end-to-end learning [40]. More generally, a GNN-based approach for graph classification is a composition of functions (or layers) of two types: 1) GNN layers (also called graph convolutional layers), which produce node representations based on the graph structure and node features (or initial node representations), and 2) pooling layers that, depending on their role in the architecture, either coarsen the network into a smaller one or produce a graph representation that can be used in a NN (also called a read-out layer) for end-to-end learning. Figure 3 provides a schematic illustration.

Here we present several GNN-based approaches suitable to use with AdaDetect, chosen for their popularity, theoretical justification and interpretability.

• GIN by [40]. In general, Graph Isomorphism Network (GIN) denotes a type of GNN layer, in which the new node representations are obtained by

  $$\begin{aligned} X' = \text {MLP} \Big ( (A+ I) X \Big ), \end{aligned}$$

  where MLP is a multilayer perceptron, A an adjacency matrix, I the identity matrix, and X is a node representation. For the task of graph classification, [40] propose to combine several GIN layers with a read-out layer that consists in summing up node representations. Moreover, the authors prove that GIN has favorable theoretical properties, namely that GIN is as powerful as the 1-WL test.

• DiffPool by [41]. GIN and most other GNNs are rather flat networks, so that they can only capture local patterns. To learn graph properties on a higher level, differentiable graph pooling (DiffPool) combine GNN layers with pooling layers that successively coarsen the graph. Coarser graphs may represent more global features of the initial graph. In each pooling operation of DiffPool, a new (coarsened) graph and new node representations are obtained by

  $$\begin{aligned} A'&= S^T A S, \quad S = \text {softmax} \Big ( \text {GNN}_{\text {pool}}(A, X) \Big ) \\ X'&= S^T \text {GNN}_{\text {embed}}(A, X) \end{aligned}$$

  where $\text {GNN}_{\text {pool}}$ and $\text {GNN}_{\text {embed}}$ are GNN layers (e.g. GIN layers). The matrix S represents a (differentiable) clustering of the nodes that is used to coarsen the input graph A. A final vector-valued graph representation is obtained using a standard read-out layer.

• DGCNN by [39]. To extend convolution neural networks (CNN) to graph-structured input, [39] introduce SortPool, a special kind of read-out layer that also acts as a coarsening operation. The SortPool layer produces a sorted representation of the nodes, such that applying classical one-dimensional CNN layers to these representations makes sense. Moreover, SortPool unifies graph sizes by truncating/extending all sorted representations to the same length, say k (where k is e.g. such that 50% of the graphs have more than k nodes). Dynamic graph CNN (DGCNN) refers to the architecture that results from combining GNN layers with SortPool and one-dimensional CNN layers.

On the one hand, GNNs have more flexibility in that they can easily take into account various characteristics of the network data, such as node attributes of any type, node correspondance, or directed edges. Moreover, graph kernels suffer from a certain computational burden in terms of the number of networks at hand. On the other hand, GNNs inherit from the lack of interpretability of NNs. Finally, while GNNs support end-to-end learning (whereas graph kernels produce fixed embeddings), their expressivity power is limited, since concerning the task of distinguishing networks GNNs cannot be more powerful than the 1-WL algorithm [40, 44].

Fig. 3

A typical GNN architecture for graph classification: N denotes the number of nodes of the input graph (A, X), L and $L'$ denote the number of GNN and pooling layers in the architecture, respectively, without the read-out layer. Each GNN layer applies a non-linear transformation on the current graph representation $A^l$ and the current node embedding matrix $X^l$ (initially A, X) and produces an updated node embedding matrix $X^{l+1}$, whereas each pooling layer coarsens the graph representation $A^l$ (initially A) into $A^{l+1}$. The final layer (read-out layer) takes the embedding matrix $X^{L+L'}$ and transforms it into a one-dimensional embedding

Results

Real complex graphs

Table 1 Summary of the used datasets

In this section, we apply the method introduced above to real datasets. For the purpose of illustration, we use graph classification datasets, where all data are labeled. A summary of the considered datasets downloaded from [47] is given in Table 1. The datasets DD and PROTEINS encode information about macromolecules. The nodes of PROTEINS represent secondary structure elements, and an edge exists if they are neighbours along the amino acid sequence. In DD, the nodes are amino acids and the edges refer to the spacial proximity. NCI1 represents chemical compounds where nodes are atoms and edges are bonds between the atoms. This dataset is relative to the cell lung cancer task. AIDS represent molecular compounds from the Antiviral Screen Databased of Active Compounds.

In our novelty detection task, we consider the graphs of one class as anomalies, and the remaining observations as references. In each task, we construct test samples and training samples by subsampling the dataset. We choose the test samples ${{\,\mathrm{\mathcal {G}_{\text {test}}}\,}}$ such that half of its networks are novelties, that is $\vert {{\,\mathrm{\mathcal {I}_{\textrm{ref}}}\,}}\vert / \vert {{\,\mathrm{\mathcal {I}_{\textrm{nov}}}\,}}\vert = 0.5$, and the size of the test samples $m=|{{\,\mathrm{\mathcal {G}_{\text {cal}}}\,}}|$ as given in Table 2. The size of the reference sample is $|{{\,\mathrm{\mathcal {G}_{\text {ref}}}\,}}|=n=2m$ and for the calibration set $\vert {{\,\mathrm{\mathcal {G}_{\text {cal}}}\,}}\vert =m$.

We apply AdaDetect with each of the graph classification methods described in the previous section: the GNN-based approaches GIN, DGCNN, and DiffPool, and one graph kernel-based approach, using the WL kernel, leading to the procedures AdaDetect-GIN, AdaDetect-DGCNN, AdaDetect-DiffPool and AdaDetect-WL. For each GNN, the architecture consists of 3 layers and 32 neurons and we train for 10 epochs with a learning rate of 0.001, and the WL kernel is used with 5 iterations. Moreover, we compare our results to the conformal anomaly detection (CAD) procedure proposed by [25], using one-class classification approaches: the one-class classifier given by the Support Vector Data Description (SVDD) method introduced in [48] with a family of functions given by either GIN, DiffPool, or DGCNN, and a one-class SVM using the WL kernel, which gives the procedures CAD-GIN, CAD-DGCNN, CAD-DiffPool, and CAD-WL.

Table 2 Settings

The FDR and TDR for the methods are evaluated on 100 subsampled data sets and the results are reported in Table 3. First, we observe that all methods control the FDR at the nominal level $\alpha =0.2$. Most often the AdaDetect version achieves a larger FDR than its CAD counterpart. Concerning the TDR, values vary a lot over the four settings and the different procedures, ranging from 0.00 to 0.95. On the data sets DD, PROTEINS and AIDS, AdaDetect with any GNN classifier outperforms the corresponding CAD version (with one exception), illustrating an important gain in power due to the AdaDetect approach. This is in line with the properties of AdaDetect reported in [27]. Concerning WL, depending on the setting, CAD is slightly doing better than AdaDetect in terms of power. The data set NC11 appears to be an inherently difficult setting as none of the procedures detects many novelties and also the FDRs are far below the nominal level $\alpha$.

Table 3 Performance of different methods on different data sets in terms of FDR (top) and TDR (bottom) with nominal level is $\alpha =0.2$. Mean values and standard deviations (in parentheses) over 100 subsampled data sets

Fig. 4

Data represented as graphs: FDR (upper row) and TDR (lower row) for AdaDetect procedures with $\alpha =0.1$ (dashed green line) for the setups described in Table . Each boxplot is based on 50 data sets. Red points indicate mean values

Fig. 5

Results of the experiments based on the extracted topological features: FDR (upper row) and TDR (lower row) for AdaDetect procedures with $\alpha =0.1$ (dashed green line) for the setups described in Table . Each boxplot is based on 50 data sets

Metabolic networks of bacteria

In this section AdaDetect is applied to a database of metabolic networks. To illustrate the performance of AdaDetect, we construct several novelty detection setups using different characteristics of the bacteria as class labels and compute the associated FDR and TDR.

Data description

The data set contains 5610 prokaryotic species from the KEGG database [4]. Among them, there are 301 archaea and 5309 bacteria. The taxonomic information was obtained from the National Center for Biotechnology Information (NCBI) Taxonomy database [49]. The reconstructed networks were provided by [8]. All information on the species was extracted in November/December 2019. The following characteristics of the bacteria are provided, which we use to build different novelty detection tasks:

• Oxygen tolerance: 917 Aerobe, 782 Facultative anaerobe, 532 Anaerobe

• Habitat: 554 Symbionts, 395 Environment, 235 Mixed

In Table 4 we provide the number of nodes, number of edges, degree assortativity, density, and average degree connectivity for the classes of bacteria. We provide the average values with the standard deviation per class. Note that the number of nodes is the same here, since we followed the network reconstruction procedure used in [7, 8] where the networks are reconstructed using the same set (the union) of nodes.

Table 4 Topological characteristics of the bacteria dataset

Experimental setup

We use the characteristics provided above to define groups of bacteria, which are then used to construct several novelty detection tasks, by labelling a bacterium as either a reference or a novelty depending on which group it belongs to. Table 5 describes which groups are considered as references or novelties in each setup.

Table 5 Description of the novelty detection tasks considered for the metabolic network dataset

In each task, we construct test samples and training samples by using the complete set of novelties and a random subset of the references as test observations, with the remaining references used for the training sample: the sample sizes $\vert {{\,\mathrm{\mathcal {I}_{\textrm{ref}}}\,}}\vert , \vert {{\,\mathrm{\mathcal {I}_{\textrm{nov}}}\,}}\vert , n$ are given in Table  for each scenario. We set $\alpha =0.1$ and use $\vert {{\,\mathrm{\mathcal {G}_{\text {cal}}}\,}}\vert = \vert {{\,\mathrm{\mathcal {G}_{\text {train}}}\,}}\vert = n/2$ for splitting the training sample.

When applying AdaDetect to these data, we observed that results are unstable and depend on the random split of the reference set ${{\,\mathrm{\mathcal {G}_{\text {ref}}}\,}}$ into subsets ${{\,\mathrm{\mathcal {G}_{\text {train}}}\,}}$ and ${{\,\mathrm{\mathcal {G}_{\text {cal}}}\,}}$. This indicates that the sample size of the reference set is small compared to the variability of the networks in the reference set. To solve this instability issue, we choose to apply each method 10 times with 10 different splits of ${{\,\mathrm{\mathcal {G}_{\text {ref}}}\,}}$ and consider the union of all detected networks as the final set of detections. We tested different number of layers, neurons, and epochs in our deep architectures. We found out that the optimal values are 3 layers with 32 neurons. We performed 10-fold cross-validation to fix these parameters. Note that we aimed to do sustainable learning, to keep the deep architecture as flat as possible, without degrading its performance, and to avoid unnecessary re-training and excessive number of epochs. We are aware of the stability and variability problems, however, in real applications, we are limited to a number of data provided by a study. Speaking about imbalanced dataset, we are completely aware of this problem. To overcame the problem of the imbalanced data, for each run of our experiments, we select the same number of observations from each class of bacteria.

Note that metabolic networks are too large for the WL algorithm, so we only consider AdaDetect with the three GNNs. The FDR and TDR are displayed in Fig. 4 for 50 randomly constructed samples. We observe that the FDR is controlled (or close to) in all settings for all methods. Concerning the power, GIN makes only very few detections, while AdaDetect with DiffPool and DGCNN are powerful procedures and depending on the setting, one or the other achieves a better TDR. We tested the state-of-the-art GNN. The relatively poor results of the GIN can be explained as follows: the GIN is unstable in tasks with big graphs [50]. The DiffPool and DGCNN are more robust to huge structures, and the explanation is probably in the pooling layers present in both architectures [51].

We also performed the experiments using the topological properties of the graphs. We extracted the topological features from the adjacency matrices, we performed training (using the cross-validation and choosing the optimal hyper-parameters) and we tested on this newly created data set. For the choice of the topological features, we got inspired by the topological features considered in [5]. Figure 5 illustrates our results using the topological features: average degree connectivity, degree centrality, betweenness centrality, degree assortativity coefficient, edge density, and average hierarchical flow. The classification methods are Random Forest (RF) with 100 estimators and depth of trees equal to 10, MLP (Multilayer Perceptron) with 3 layers and 5 neurons in each layer, and Gradient Boosting (GradBoost) with 100 estimators. We observed that the deep learning architectures, DiffPool and DGCNN reach a better performance. We conclude that taking into account the whole graph is beneficial in the current application, since the resulting estimated probabilities to belong to a class are estimated more accurately compared to a dataset based on the topological features.

Metabolic networks of cancer patients

We test the proposed approach and run several experiments on the rich recent dataset TumorMet [19]. In particular, we applied our method on two tasks of the repository, namely, on Lung data (TCGA-LUAD, 428 patients; TCGA-LUSC, 400 patients) and Kidney data (TCGA-KIRC, 483 patients; TCGA-KIRP, 253 patients). We focused on the metabolites-based networks where each network corresponds to a patient, the edges are weighted using the gene expression. Therefore, the graphs are directed and weighted. More detailed description of the data and the topological properties of the graphs can be found in [19]. Figure 6 illustrates the results of the numerical experiments on the Lung and Kidney tasks. We can conclude that although our results are positive, they also confirm that controlling FDR and maximising TDR are challenging problems, whose results heavily rely on a classification model and on the estimated probabilities to belong to classes.

Fig. 6

Results of the numerical experiments on Lung and Kidney tasks from the repository of tumor metabolic networks [19]

Discussion and Conclusions

Conformal prediction is an emerging direction in medical and biological applications, since statistical machine learning models have to be developed and applied with caution, especially in high-stake domains.

We propose a powerful tool, applicable to complex structures, to control the desired risk level. To our best knowledge, we are the first ones to challenge this task. Although our method achieves remarkable performance, there is still room for further research. So, our next step is to increase the interpretability of graph embeddings used in the method. Note that different deep or flat models lead to various graph embeddings or representations, and it is an important question how to find an optimal embedding. Particularly, the metabolic networks are huge complex graphs, and their reconstruction and representation can vary according to the scientific aim, e.g., some ubiquitous metabolites can be omitted, enzymes can be included or not, an algorithm of network reconstruction can easily result in different graphs; the direction of reactions is not unique and can also vary according to an ontology used for the graph reconstruction.

The unique role of our approach is the control of the desired risk level of prediction, i.e., our method returns results with guarantees on the classification quality. Such a guarantee is absolutely needed to control uncertainty in the prediction models: if a model is uncertain, it is better to abstain from taking decision. The task to cope with the uncertainty is even more challenging for complex data, such as metabolic networks. In our contribution, we show that our method controls the empirical error at a desired (5%, 10%, or 20%, depending on a task) level.

Note that our conformal prediction approach can be used in combination with any classifier. The reason to consider GNNs here is that graph features are not chosen manually, but graph representations are learned automatically from the data by the GNN.

An open ambitious question is how to relate the data representation and the FDR control, and whether a unified efficient framework can be proposed and developed. Another general question is the anomaly detection with techniques different from the conformal prediction. So, in recent years the detection of adversarial examples was identified as an important issue, and [52] proposed to tackle this problem with ensemble approaches. This setting is similar to ours in some aspects, however, the ensemble method of [52] makes its final decision based on a vector of scores. Note that our method is a statistically sound method based on a non-conformity score. One of the promising research avenues is to study how to combine several scores (e.g., from an ensemble method). The answer to this problem is not straightforward, since it is not clear how to get statistical guarantees with an ensemble approach.

We also plan to apply the method to different tasks and, in particular, to medical challenges. We started to explore the recent rich repository TumorMet [19] dedicated to various pathologies and that contains different metabolic graph representations. Such datasets have to be explored more thoroughly.

Finally, in this contribution is focused on the confident classification. However, much more deeper analysis of the identified novelties in graphs and tight collaboration with human experts in a specific domain is a further step to understand mechanisms of phenomena, e.g., human pathologies.