11institutetext:  Department of Electronics and Computer Science
Pontificia Universidad Javeriana, Cali, Colombia 11email:
<EMAIL_ADDRESS>
# Feature extraction with Spectral Clustering for Gene Function Prediction
using Hierarchical Multi-label Classification

Miguel Romero  Oscar Ramírez  Jorge Finke  Camilo Rocha

###### Abstract

Gene annotation addresses the problem of predicting unknown associations
between gene and functions (e.g., biological processes) of a specific
organism. Despite recent advances, the cost and time demanded by annotation
procedures that rely largely on in vivo biological experiments remain
prohibitively high. This paper presents a novel in silico approach for to the
annotation problem that combines cluster analysis and hierarchical multi-label
classification (HMC). The approach uses spectral clustering to extract new
features from the gene co-expression network (GCN) and enrich the prediction
task. HMC is used to build multiple estimators that consider the hierarchical
structure of gene functions. The proposed approach is applied to a case study
on Zea mays, one of the most dominant and productive crops in the world. The
results illustrate how in silico approaches are key to reduce the time and
costs of gene annotation. More specifically, they highlight the importance of:
(i) building new features that represent the structure of gene relationships
in GCNs to annotate genes; and (ii) taking into account the structure of
biological processes to obtain consistent predictions.

##  1 Introduction

Identifying the association of genes to functions is key to gain insight into
how genomes serve as blueprints for life, e.g., to develop treatments for
specific conditions or enhance tolerance to environmental stresses [29, 34,
37]. Numerous studies have used co-expression data to predict specific
biological functions and processes [23, 26, 31, 33]. Intuitively, genes are
reported to co-express whenever they are simultaneously active, which suggests
that they are associated to one or more common biological processes.

Under this hypothesis, characterizing gene interactions as a gene co-
expression network (GCN) may assist to identify unknown functional annotations
in a genome. Co-expression networks are generally represented as undirected
weighted graphs, where vertices denote genes and weighted edges indicate the
strength of the co-expression between two genes. A detailed analysis of the
structure and distribution of gene relationships in GCNs provides additional
clues that facilitate the prediction of gene functions [32].

However, the cost and time requirements to annotate genes using in vivo
biological experimentation remains prohibitively high [4, 42]. To overcome
this limitation, hybrid approaches that integrate existing knowledge of gene-
function associations and in silico methods have been proposed [5, 7, 17, 27].
While they have shown great promise, given the extreme combinatorial nature of
the problem, annotating genes in an efficient manner remains an open
challenge.

Functional annotations are defined by the Gene Ontology (GO), which contains
three main types of annotations: biological processes, molecular functions,
and cellular component [8]. These annotations, commonly known as GO terms, are
structured in a hierarchy and defined as a directed acyclic graph (DAG). Gene
annotation approaches generally ignore the relationships among biological
processes, even though these relationships are key to improve the accuracy and
avoid inconsistency in predictions. A prediction is said to be inconsistent
w.r.t. the GO hierarchy when a gene is inferred to have a particular function
$a$, but it is not inferred to have all ancestor of $a$. In other words, an
inconsistent prediction states that the prediction does not satisfy the
ancestral relations between GO terms. Satisfying ancestral constraints is
often referred to as the true-path rule in GO [32, 1] and as the hierarchical
constraint in HMC [35].

This paper presents a feature extraction approach for in silico annotation of
genes. It follows a network-based approximation that uses cluster analysis and
hierarchical multi-label classification (HMC) for building a predictor that
assigns functions to genes satisfying the true-path rule. Cluster analysis
plays the role of enriching the information available for predicting gene-
function associations by extracting new features that represent structural
properties of the GCN. That is, co-expression relations are used to identify
gene clusters that ultimately help in associating functions to genes (i.e.,
guilt by association, see [24]). It has been shown in [28] that new features
built from the GCN and associations between genes and functions with the
spectral clustering algorithm are key to improve the prediction performance in
the gene annotation problem. The results in [28] show that using other
features associated to structural properties of the GCN and gene functional
information lead to lower performance.

Furthermore, the extracted features are filtered (using SHAP) based on their
impact in the prediction task and HMC is used to predict gene-function
associations that take into account the relations between biological
functions. The proposed approach illustrates how the performance of gene
annotation is improved by combining: (i) new information extracted from the
GCN; and (ii) classification methods that consider the relation between gene
functions.

This approach is applied to a case study on Zea mays, one of the most dominant
and productive crops. Zea mays serves a variety of purposes, including animal
feed and derivatives for human consumption and ethanol [41]. The co-expression
information used in the study is imported from the ATTED-II database [21]. The
resulting GCN, modeled as a weighted graph, comprises 26 131 vertices (i.e.,
genes) and 44 621 533 edges. The functional information (i.e., known gene-
function associations) is taken from DAVID Bioinformatics Resources [10]. It
contains a total of 255 865 annotations of biological processes for maize,
i.e., pathways to which a gene contributes. The results highlight the
importance of extracted features that represent structural properties of the
GCN and the hierarchical structure of biological processes with HMC to improve
prediction performance. Ultimately, the results provide experimental (in
silico) evidence that the proposed approach is a viable and promising
approximation to gene function prediction.

This paper is a significant extended version of [28] that:

  * •

Addresses the gene function prediction as a hierarchical multi-label
classification problem by considering the structure of gene functions. That is
ancestral relationships are represented as a DAG [8].

  * •

Analyzes a larger functional database for the case study of maize. The number
of genes associated to at least one function increased from 5 361 to 10 049.
The new dataset consists of 255 865 associations between genes and functions,
and 7 021 relations between functions.

  * •

Concludes that the ancestral relations between functions and the features
extracted from the GCN improve the prediction performance in the gene function
prediction task when addressed as a hierarchical multi-label classification
problem.

The remainder of the paper is organized as follows. Section 2 reviews some
preliminaries. Section 3 introduces the approach to extract features from the
gene co-expression network using cluster analysis. The proposed approach to
predict gene functions, based on hierarchical multi-label classification is
presented in Section 4. Section 5 presents the case study for the Zea mays
species. Finally, Section 6 draws some concluding remarks and future research
directions.

##  2 Preliminaries

This section presents preliminaries on spectral clustering, gene co-expression
networks, gene function prediction, hierarchical multi-label classification,
and SHAP feature contribution.

###  2.1 Spectral clustering

The aim of applying cluster analysis on a network is to identify groups of
vertices sharing a (parametric) notion of similarity [39, 25]. Usually,
distance or centrality metrics are used for clustering. Spectral clustering is
a clustering method with foundations in algebraic graph theory [11]. It has
been shown that spectral clustering has better overall performance across
different areas of applications [20]. Given a graph $G$, the spectral
clustering decomposition of $G$ can be represented by the equation
$\mathbf{L}=\mathbf{D}-\mathbf{A}$, where $\mathbf{L}$ is the Laplacian,
$\mathbf{D}$ is the degree (i.e., a diagonal matrix with the number of edges
incident to each node), and $\mathbf{A}$ the adjacency matrices of $G$.
Spectral clustering uses, say, the $n$ eigenvectors associated to the $n$
smallest nonzero eigenvalues of $\mathbf{L}$. In this way, each node of the
graph gets a coordinate in $\mathbb{R}^{n}$. The resulting collection of
eigenvectors serve as input to a clustering algorithm (e.g., k-means) that
groups the nodes in $n$ clusters.

###  2.2 Gene Co-expression Network

A gene co-expression network (GCN) is represented as an undirected graph where
each vertex represents a gene and each edge the level of co-expression between
two genes.

###### Definition 1

Let $V$ be a set of genes, $E$ a set of edges that connect pairs of genes, and
$w:E\to\mathbb{R}_{\geq 0}$ a weight function. A _(weighted) gene co-
expression network_ is a weighted graph $G=(V,E,w)$.

The set of genes $V$ in a co-expression network is particular to the genome
under study. The correlation of expression profiles between each pair of genes
is measured, commonly, using the Pearson correlation coefficient. Every pair
of genes is assigned and ranked according to a relationship measure, and a
threshold is used as a cut-off value to determine $E$. The weight function $w$
denotes the strength of the co-expression between each pair of genes in $V$.
For example, in the ATTED-II database, the co-expression relation between any
pair of genes is measured as a $z$-score expressed as a function of the co-
expression index LS (Logit Score) [21, 22].

###  2.3 Gene Function Prediction

In an annotated gene co-expression network, each gene is associated with the
collection of biological functions to which it is related (e.g., through in
vivo experiments).

###### Definition 2

Let $A$ be a set of biological functions. An _annotated gene co-expression
network_ is a gene co-expression network $G=(V,E,w)$ complemented with an
annotation function $\phi:V\to 2^{A}$.

The problem of predicting gene functions can be explained as follows. Given an
annotated co-expression network $G=(V,E,w)$ with annotation function $\phi$,
the goal is to use the information represented by $\phi$, together with
additional information (e.g., features of $G$), to obtain a function
$\psi:V\to 2^{A}$ that extends $\phi$. Associations between genes and
functions not present in $\phi$ have either not been found through in vivo
experiments, or do not exist in a biological sense. The new associations
identified by $\psi$ are a suggestion of functions that need to be verified
through in vivo experiments. The function $\psi$ can be built from a predictor
of gene functions, e.g., based on a supervised machine learning model.

###  2.4 Hierarchical Multi-label Classification

Node classification refers to the task of predicting a node class for an input
data based on the information of other nodes in the network [2]. In general,
node classification problems can be categorize into three different types:
binary classification refers to predict one attribute (target) with two
classes (for example, positive and negative) [12]; multi-class classification
refers to the case where the attribute to be predicted has more than two
classes and are mutually exclusive (for example, the brand of a car) [18]; and
multi-label classification refers to predicting an attribute with at least two
classes, but where an instance could be associated to more than one class (for
example, the gene function prediction problem) [36].

Although the aforementioned prediction methods are frequently used, they do
not consider hierarchical relations between classes. For such scenarios,
hierarchical multi-label classification (HMC) addresses the task of structured
output prediction where the classes are organized into a hierarchy and an
instance may belong to multiple classes. In many problems, such as gene
function prediction, classes inherently satisfy these conditions [14]. Authors
in [30] expose that there are two types of methods to explore the hierarchical
structure. First, top-down or local classifiers refer to partially predict the
classes in the hierarchy from the top to the bottom. Second, big-bang or
global classifiers refer to use a single classifier that considers the entire
hierarchy at once.

Classifiers that ignore the class relationships, by predicting only the leaf
classes in the hierarchy or predicting each class independently, often lead to
_inconsistent predictions_. This refers to the fact that a node is inferred to
have a particular class $a$, but the outcome of the classifier fails to infer
the node’s association to all ancestor classes of $a$ in the hierarchy. In
other words, an inconsistent prediction states that the prediction does not
satisfy the hierarchy for some class $a$. Satisfying ancestral constraints is
often referred to as the true-path rule in GO [32, 1] and as the hierarchical
constraint in HMC [35].

Figure 1 illustrates the four HMC methods used in this work: Local classifier
per node (lcn) consists of training one binary classifier for each class in
the hierarchy except the root. Local classifier per parent node (lcpn)
consists of training a multi-label classifier for each parent node in the
hierarchy to distinguish between its child classes. Local classifier per level
(lcl) consists of training one multi-label classifier for each level of the
class hierarchy except for the root. Global classifier consists of building a
single multi-label classifier taking into account the hierarchy as a whole
during a single run. The global classifier can assign classes at potentially
every level of the hierarchy to an instance.

Figure 1: Example of global and local methods for hierarchical multi-label
classification. Given a hierarchy of classes ($r$, $a$, $b$, $c$, $d$, $e$,
and $f$), the dashed boxes show the number of classifiers required for each
method. Note that the lcn, lcpn, lcl, and global classifiers require 6, 4, 3,
and 1 predictors, respectively.

###  2.5 SHAP feature contribution

The performance of classification algorithms is partly determined by the
features used to train a particular predictor. SHAP (SHapley Additive
exPlanation) is a framework that computes the importance values for each
feature in a dataset using concepts from game theory [15, 16]. SHAP assigns
Shapely values to explain which features in the model are the most important
for prediction by calculating the changes in the prediction when features are
conditioned. Given a predictor and a training set, SHAP computes a matrix with
the same dimensions of the predictor’s output containing the Shapely values
for each instance and class. For example, in a binary classification problem
and a training set of $n$ instances, the output of SHAP is a matrix of
dimension $n\times 2$ (there are two classes, positive and negative). In
multi-label classification problems, the output is a matrix of dimension
$n\times 2$ for each class, since classes are not mutually exclusive and the
outcome is either positive or negative for each class.

##  3 Clustering-based Feature Extraction

The approach for extracting features from the GCN using a clustering algorithm
and Gene Ontology term enrichment is presented. It combines information from
the GCN, and the associations between genes and functions to create features
capturing topological properties of the GCN.

The inputs of the approach are a GCN, denoted by $G=(V,E,w)$, a set of
(biological) functions $A$, an annotation function $\phi:V\to 2^{A}$, and a
set $K=\\{k_{0},\dots,k_{m-1}\\}$ for sampling the number of clusters. The
annotation function $\phi$ must satisfy true-path rule for the GO hierarchy
[1, 32]. That is, if a gene is associated to a function, then it must also be
associated to every ancestor of the function in the hierarchy, and if a gene
is not associated to a function, then it must not be associated to any of its
descendants.

The outputs are two feature matrices $J_{G}$ and $J_{F}$, of dimension
$V\times A\cdot K\rightarrow[0,1]$, specifying the likelihood of the genes $V$
to be associated to the functions in $A$ when the graph is decomposed in $m$
clusters. Matrices $J_{G}$ and $J_{F}$ correspond to the GCN (that is the
graph $G$) and an affinity graph defined the next subsection.

The feature extraction approach consists of three stages, which are depicted
in Figure 2. First, an affinity graph $F$ with information in $\phi$ is
created from $G$. Second, the spectral clustering algorithm is applied to both
$G$ and its enriched version $F$ for the $m$ different number of clusters
specified in $K$. Third, the Gene Ontology term enrichment technique is used
to create $m$ features for each function $a\in A$, corresponding to the number
of clusters in $K$.

Figure 2: The clustering-based feature extraction approach consists of three
stages. Namely, creation of affinity graph, clustering computation, and Gene
Ontology term enrichment. Its inputs are a GCN, denoted by $G=(V,E,w)$, a set
of functions $A$, an annotation function $\phi:V\to 2^{A}$, and a set
$K=\\{k_{0},\dots,k_{m-1}\\}$. Its output are two feature matrices (for both
$G$ and its enriched version $F$) of dimension $V\times A\cdot
K\rightarrow[0,1]$ that specify how likely it is for the genes to be
associated to the functions in $A$ when the graph is decomposed $m$ clusters,
each of size $k_{i}$, for $0\leq i\leq m$.

###  3.1 Affinity Graph Creation

An affinity graph $F=(V,E,w_{F})$ between $G$ and $\phi$ is built. Its weight
function is defined as the mean between the co-expression weight specified by
$w$ and the proportion of shared functions between genes specified by $\phi$.

###### Definition 3

The _weight function_ $w_{F}:V\times V\to[0,1]$ is defined for any $u,v\in V$
as

$\displaystyle
w_{F}(u,v)=\frac{1}{2}\left(\frac{w(u,v)-1}{\max(w)-1}+\frac{|\phi(u)\cup\phi(v)|}{|\phi(u)\cap\phi(v)|}\right),$

where $\max(w)$ denotes the maximum value in the range of $w$ (which exists
because $w$ is finite).

Under the assumption that at least one element in the range of $w$ is greater
than 1, it is guaranteed that the range of $w_{F}$ is $[0,1]$ (because
$w:V\times V\to[1,\infty)$). This is indeed the case, in practice, because the
co-expression between two genes in the GCN is quantified in terms of the
$z$-score, which is highly unlikely to be 1 for all pairs of genes.

###  3.2 Gene Clustering

The spectral clustering algorithm is applied independently to each graph
$X\in\\{G,F\\}$ to decompose $X$ (i.e., group the genes $V$) using the number
of clusters specified by $K=\\{k_{0},\dots,k_{m-1}\\}$. The decomposition of
$X$ is performed $m$ times, once per $k$ in $K$. The adjacency matrices of the
weighted and undirected graphs $G$ and $F$ are used as the precomputed
affinity matrices required for the spectral clustering algorithm. The outcome
of the clustering algorithm is an assignment from nodes to clusters of size
$k$, for each $k\in K$. More precisely, the outputs of this stage are the
matrices $I_{X}:V\times K\to[0,1]$, where each column $0\leq i<m$ represents
the decomposition of $X$ in $k_{i}$ clusters.

###  3.3 Gene Enrichment

The goal of this stage is to produce a matrix $J_{X}:V\times A\cdot K\to[0,1]$
for each $X\in\\{G,F\\}$, specifying how likely it is for the genes to be
associated to every function $a\in A$ when $X$ is decomposed in the given
number of clusters.

For each decomposition from the previous stage (i.e., each column of the
matrices $I_{X}$) and function $a\in A$, the resulting clusters are used to
compute whether a significant number of members associated to function $a$ is
(locally) present. Intuitively, if genes that are grouped together have a
strong co-expression relation and most of the group are associated to gene
function $a$, then the remaining genes are also likely to be associated to $a$
(i.e., guilt by association, see [24]). In this way, for each $v\in V$, $a\in
A$, and $k\in K$, the entry $J_{X}(v,a\cdot k)$ is a $p$-value indicating if
the function $a$ is over-represented in the decomposition of $k$ clusters of
$X$. This process is commonly known as Gene Ontology term enrichment and may
use different statistical tests, such as, Fisher’s exact test [38].

##  4 Hierarchical Multi-label Classification for Gene Function Prediction

This section presents the approach for gene function prediction using HMC to
create a predictor, enriched with the information of the features created in
Section 3.

The GO hierarchy is defined as a directed acyclic graph (DAG) containing three
main types of annotations: biological processes, molecular functions, and
cellular component [8]. This work focuses on biological processes, i.e., a
subgraph of the GO hierarchy that contains 28 roots (i.e., functions in the GO
hierarchy with null indegree). This subgraph is denoted as $H=(A,R)$, where
$A$ is the set of biological processes and $R$ the binary relation
representing ancestral relations between pairs of biological processes (i.e.,
$(a,b)\in R$ means that function $b$ is ancestor of function $a$ in the GO
hierarchy). The topological-sorting traversal algorithm presented in [27] is
used to transform the GO hierarchy of biological processes into a tree. As a
result, the hierarchy is split into several components, i.e., subtrees of $H$
called sub-hierarchies. Each sub-hierarchy,
$H^{\prime}=(A^{\prime},R^{\prime})$ with $A^{\prime}\subseteq A$,
$R^{\prime}\subseteq R$, and $r\in A^{\prime}$ the root, is associated to a
subgraph $G^{\prime}=(V^{\prime},E^{\prime},w)$ containing all genes $v\in V$
associated to $r$, i.e., $V^{\prime}=\phi^{-1}(r)$. Note that, the proposed
approach is independently applied to each sub-hierarchy.

The inputs of the approach are a sub-hierarchy
$H^{\prime}=(A^{\prime},R^{\prime})$, a subgraph of the GCN, denoted by
$G^{\prime}=(V^{\prime},E^{\prime},w)$, where $V^{\prime}\subseteq V$ and
$E^{\prime}\subseteq E$, an annotation function $\phi:V\to 2^{A^{\prime}}$,
the matrices $J_{G}$ and $J_{F}$ resulting from Section 3, and a constant
value $c\in[0,1]$ for feature selection. The output is a function
$\psi:V^{\prime}\times A^{\prime}\to[0,1]$, specifying, for each gene $v\in
V^{\prime}$, the probability $\psi(v,a)$ of $v$ being associated to function
$a\in A^{\prime}$.

First, sub-matrices $J_{G}^{\prime}$ and $J_{F}^{\prime}$ are created from
$J_{G}$ and $J_{F}$, by respectively considering only the genes
$V^{\prime}\subseteq V$ and functions $A^{\prime}\subseteq A$. These sub-
matrices represent structural properties of the GCN subgraph $G^{\prime}$, and
associations between genes and functions based on multiple partitions of each
graph. Figure 3 illustrates the prediction approach. The reminder of this
section is devoted to detailing the prediction approach.

Figure 3: The prediction approach mainly consists of two stages, feature
selection with SHAP and hierarchical multi-label classification. Its inputs
are a sub-hierarchy $H^{\prime}=(A^{\prime},R^{\prime})$, a subgraph of the
GCN $G^{\prime}=(V^{\prime},E^{\prime},w)$, an annotation function $\phi:V\to
2^{A^{\prime}}$ that satisfy the sub-hierarchy $H^{\prime}$, the sub-matrices
of $J_{G}$ and $J_{F}$ containing only the functions $A^{\prime}$ and genes
$V^{\prime}$, and a constant value $c\in[0,1]$ for feature selection. Its
output is a function $\psi:V^{\prime}\times A^{\prime}\to[0,1]$, which
indicates for each gene $v\in V^{\prime}$, the probabilities $\psi(v,a)$ of
$v$ being associated to function $a\in A^{\prime}$.

SHAP filters the extracted features with more impact in the prediction task,
and HMC is used to predict associations between genes and functions without
inconsistencies (i.e., complying the true-path rule). Since local HMC methods
use more than one predictor per hierarchy, the feature selection is executed
for each predictor independently, considering only the features related to the
functions being predicted, denoted by $A^{\prime\prime}\subseteq A^{\prime}$.
For example, consider the function hierarchy and a local classifier per level
method depicted in Figure 4. The predictor for level 1 predicts functions $a$
and $b$, so only the features associated to functions $a$ and $b$ are
considered for the feature selection.

Figure 4: Gene function prediction considering the function hierarchy and
using a local classifier per level method. The predictor for the level 1
predicts functions $a$ and $b$, so only the features from $J_{G}^{\prime}$ and
$J_{F}^{\prime}$ associated to functions $a$ and $b$ are considered for the
feature selection.

###  4.1 Feature Selection

The aim of feature selection is to produce a matrix
$J:V^{\prime}\times\Theta(c)\to[0,1]$ by selecting a reduced number of
significant features from $J_{G}^{\prime}$ and $J_{F}^{\prime}$. The number of
selected features is denoted by $0\leq\Theta(c)\leq
2m\cdot|A^{\prime\prime}|$, where $m\cdot|A^{\prime\prime}|$ is the number of
features in each matrix $J_{G}^{\prime}$ and $J_{F}^{\prime}$, denoted as $q$
(that is $q=m\cdot|A^{\prime\prime}|$).

Feature selection is conveyed from $J_{G}^{\prime}$ and $J_{F}^{\prime}$ to
$J$ using SHAP. Let $J_{G+F}^{\prime}$ denote the matrix resulting from
extending $J_{G}^{\prime}$ with the $q$ features of $J_{F}^{\prime}$. That is,
for each $v\in V^{\prime}$, the expression $J_{G+F}^{\prime}(v,\\_)$ denotes a
function with domain $[0,2q)$ and range $[0,1]$, where the values in $[0,q)$
denote the $p$-values associated to $v$ in $G$ and the values in $[q,2q)$ the
ones associated to $v$ in the enriched version of $G$. For each entry
$J_{G+F}^{\prime}(v,j)$, with $v\in V^{\prime}$ and $0\leq j<2q$, the mean
absolute SHAP value $s_{(v,j)}$ is computed after a large enough number of
Shapely values are computed (executions of SHAP). Features are selected based
on the cutoff

$c\cdot\sum_{j=0}^{2q-1}s_{(v,j)},$

i.e., on the sum of mean absolute values by a factor of the input constant
$c$. The first $\Theta(c)$ features, sorted from greater to lower mean
absolute SHAP value, are selected as to reach the given cutoff.

Note that the input constant $c$ is key for selecting the number of
significant features. The idea is to set $c$ so as to find a balance between
prediction efficiency and the computational cost of building the predictor.

###  4.2 Training and Prediction

This stage comprises a process that combines two supervised machine learning
techniques/tools to build the predictor $\psi$. In particular, stratified
$k$-fold cross-validation and hierarchical multi-label classification are used
sequentially in a pipeline.

The pipeline takes as input the matrix $J$, which specifies the significant
features of $J_{G}^{\prime}$ and $J_{F}^{\prime}$, the sub-hierarchy
$H^{\prime}$ and the annotation function $\phi$. First, $k$-fold is applied to
split the dataset into $k$ different folds for cross validation (note that $k$
is not related to the input $K$). That is, each fold is used as a test set,
while the remaining $k-1$ folds are used for training. Recall that $k$-fold
cross validation aims to overcome overfitting in training. Furthermore, one or
multiple random forest classifiers are build and used for prediction, the
number of classifiers depends on the HMC method. Randoms forest is selected
for this approach since it is a tree-based and multi-label classification
algorithm, which is interpretable (SHAP can be applied). The parameter values
used for random forest classifiers, differently from the default scikit-learn
values, are: 200 estimators (n_estimators) and minimum number of samples of 5
(min_samples_split).

Additionally, some HMC methods require an extra step to keep prediction
consistent w.r.t. the sub-hierarchy $H^{\prime}$ (i.e., comply the true-path
rule). The probability of association between a function $v\in V^{\prime}$ and
a function $a\in A^{\prime}$ must be lower than the probability of association
between the same gene and the ancestor of $a$ in $H^{\prime}$. To satisfy this
constraint cumulative probabilities are computed throughout the paths in
$H^{\prime}$. That is, for each gene $v\in V$ and functions $(a,b)\in R$, the
predicted probability of the association between $v$ and $a$ is multiplied by
the predicted probability of association between $v$ and $b$ (its ancestor).
This process is repeated for every path in the hierarchy from the root to the
leaves.

The output of this stage is the predictor $\psi$, i.e., the probabilities of
associations between the genes in $V^{\prime}$ and functions $A^{\prime}$.
Note that the predictor $\psi$ satisfies the true-path rule.

###  4.3 Performance Evaluation

It is often the case in HMC datasets that individual classes have few positive
instances. In genome annotation, typically only a few genes are associated to
specific functions. This implies that for most classes (deeper in the
hierarchy), the number of negative instances by far exceeds the number of
positive instances. Hence, the real focus is recognizing the positive
instances (predict associations between genes and functions), rather than
correctly predicting the negative ones (predict that a function is not
associated to a given gene). Although ROC curves are better known, their area
under the curve is higher if a model correctly predicts negative instances,
which is not suitable for HMC problems.

For this reasons, the measures (based on the precision-recall (PR) curve)
introduced by [35] are used for evaluation.

#### Area under the average PR curve.

The first metric transforms the multi-label problem into a binary one by
computing the precision and recall for all functions $A^{\prime}$ together.
This corresponds to micro-averaging the precision and recall.

The output of the prediction stage are the probabilities of associations
between genes $V^{\prime}$ and functions $A^{\prime}$. Thereby, instead of
selecting a single threshold to compute precision and recall, multiple
thresholds are used to create a PR curve. In the PR curve each point represent
the precision and recall for a give threshold that can be computed as:

$\overline{\textnormal{Prec}}=\frac{\sum_{i}TP_{i}}{\sum_{i}TP_{i}+\sum_{i}FP_{i}},\quad\textnormal{and}\quad\overline{\textnormal{Rec}}=\frac{\sum_{i}TP_{i}}{\sum_{i}TP_{i}+\sum_{i}FN_{i}}.$

Note that $i$ ranges over all functions $A^{\prime}$, i.e., precision and
recall are computed for all functions together. The area under this curve is
denoted as AU($\overline{\textnormal{PRC}}$).

#### Average area under the PR curves.

The second metric corresponds to the (weighted) average of the areas under the
PR curves for all functions $A^{\prime}$. This metric, referred as macro-
average of precision and recall, can be computed as follows:

$\overline{\textnormal{AUPRC}}_{w_{1},w_{2},\dots,w_{|A^{\prime}|}}=\sum_{i}w_{i}\cdot\textnormal{AUPRC}_{i}.$

If the weights of all functions are the same (i.e., $1/|A^{\prime}|$) the
metric is denoted as $\overline{\textnormal{AUPRC}}$. In addition, weights can
also be defined based on the number of genes associated to functions in
$\phi$, i.e., $w_{a}=|\phi^{-1}(a)|/\sum_{i}|\phi^{-1}(i)|$ for $a\in A$. In
the later case, denoted as $\overline{\textnormal{AUPRC}_{w}}$, more frequent
functions get higher weight. Note that one point in the weighted PR curve
corresponds to the (weighted) average of the AUPRC of all functions
$A^{\prime}$ given a threshold.

##  5 Case study: Zea mays

Next section describes a case study on applying the feature extraction and
prediction approach presented in Sections 3 and 4 to maize (Zea mays). First,
the maize data used for the case study is described. Second, the proposed
approach is applied to the maize data. Lastly, the performance of the proposed
approach is compared to two models trained using each set of features $J_{G}$
and $J_{F}$, independently.

###  5.1 Data Description and Feature Extraction

The co-expression information used in the study is imported from the ATTED-II
database [21]. The gene co-expression network $G=(V,E,w)$ comprises 26 131
vertices (genes) and 44 621 533 edges. In this case, a $z$-score threshold of
$1$ is used as the cut-off measure for $G$, i.e., $E$ contains edges $e$ that
satisfy $w(e)\geq 1$ (most of them satisfying $w(e)>1$). Note that the highest
value is assigned to the strongest connections. The functional information for
this network is taken from DAVID Bioinformatics Resources [10] (2021 update);
it contains annotations of biological processes, i.e., pathways to which a
gene contributes. It is important to note that genes may be associated to
several biological processes, and biological processes may be associated to
multiple genes. The database comprises 3 924 biological processes $A$ and 7
021 ancestral relations $R$ between these functions, that represent the
hierarchy $H=(A,R)$ of the GO [8]. A total of 255 865 association between
genes and functions are considered, these associations represent the
annotation function $\phi:V\rightarrow 2^{A}$.

The feature extraction approach is applied with the inputs $G$, $A$, $\phi$
and $K=\\{10,20,\dots,100\\}$ (values are incremented in steps of 10 up to
100). The outputs are the feature matrices $J_{G}$ and $J_{F}$ that specify
how likely it is for the maize genes $V$ to be associated to the biological
processes $A$ when the graph is decomposed in the number of clusters in $K$.

Moreover, only functions associated to more than 200 genes have been
considered, so the number of functions in the resulting sub-hierarchies is
tractable regarding the dimension of the output of SHAP (see Section 2).
Recall that the Gene Ontology hierarchy splits into 28 sub-hierarchies when
considering only biological processes. Additionally, all sub-hierarchies with
less than 10 functions are discarded and the topological-sorting algorithm
introduced in [27] is used to transform the sub-hierarchies, represented as
DAGs, into trees. For each ancestral relation $(a,b)\in R$ ($b$ is ancestor of
$a$), the algorithm assigns a weight as the ratio of the number of genes
associated to the $a$ to the number of genes associated to $b$. Then, for each
function $a\in A^{\prime}$ with more than one parent, only the one with the
higher weight remains (ties are broken arbitrarily).

Root |  Description | Functions | Genes |  Functions per level
---|---|---|---|---
GO:0050896 |  response to stimulus | 13 | 1733 |  5/5/2
GO:0051179 |  localization | 25 | 1497 |  3/5/9/6/1
GO:0065007 |  biological regulation | 37 | 2647 |  2/5/11/10/4/2/2
GO:0008152 |  metabolic process | 92 | 6596 |  8/18/38/12/7/6/2
GO:0009987 |  cellular process | 92 | 8005 |  13/19/19/17/13/8/2
Table 1: Resulting sub-hierarchies $H^{\prime}$ of biological processes for
maize. The identifier and description of each root function $r$ is presented
in the first and second columns, respectively. The third column shows the
number of functions $A^{\prime}$ within each sub-hierarchy and the fourth
column shows the number maize genes in the GCN subgraph $G^{\prime}$
associated to $H^{\prime}$. The last column shows the number of functions per
level, e.g., the first sub-hierarchy has 3 levels and there are 5, 5, and 2
functions on each level.

As result, there are 5 sub-hierarchies of biological processes. Table 1
describes each sub-hierarchy $H^{\prime}$, starting by the root term $r$ and
its description, following the number of functions $A^{\prime}$ and the number
of genes $V^{\prime}$ in the associated GCN subgraph $G^{\prime}$. The
prediction approach is applied to each sub-hierarchy $H^{\prime}$
independently. The remaining input parameter for the prediction approach is
$c=0.9$ (recall that this parameter is used to filter the most relevant
features according to their mean SHAP value). Figure 5 depicts the number of
classifiers trained per HMC method and sub-hierarchy. Note that the global
method requires one classifier per hierarchy, while the lcn requires
$|A^{\prime}|-1$ classifiers.

Figure 5: number of classifiers trained per HMC method and sub-hierarchy. The
lcn requires $|A^{\prime}|-1$ classifiers. The lcpn requires as many
classifiers as functions with children in $H^{\prime}$. The lcl requires as
many classifiers as the number of levels in $H^{\prime}$. At last, the global
method requires one classifier per hierarchy.

###  5.2 Summary of Results

Figure 6: Prediction performance of the proposed approach measured with the
area under the average PR curve, i.e., AU($\overline{\textnormal{PRC}}$). The
performance is measured independently per sub-hierarchy.

Figure 6 presents the prediction performance of the proposed approach measured
with the AU($\overline{\textnormal{PRC}}$) (denoted as micro) for four HMC
methods, namely, local classifier per node (lcn), local classifier per parent
node (lcpn), local classifier per level (lcl), and global classifier. In
general, it can be seen that all methods get a high area under the average PR
curve, but the global classifier outperforms the local methods for all sub-
hierarchies. The proposed approach identifies the associations between genes
and functions by using the features extracted from the GCN $G$ and the
affinity graph $F$, and considering the ancestral relations of the biological
processes. The global method obtains the best performance, followed by the
lcpn and the lcl. Using multi-label classifiers is better than using a binary
classifier for each function, i.e., lcn method.

Figure 7: Prediction performance of the proposed approach measured with the
average area under the PR curve, i.e., $\overline{\textnormal{AUPRC}}$. The
performance is measured independently per sub-hierarchy.

The micro score measures the overall performance of all functions within a
sub-hierarchy without distinguishing between them. Figure 7 presents the
prediction performance measured with the $\overline{\textnormal{AUPRC}}$,
denoted as macro. The macro score measure the prediction performance for each
function individually and then takes the average. The conclusion is similar,
the global method outperforms the local ones.

Figure 8: Prediction performance of the proposed approach measured with the
average area under the PR curve, i.e., $\overline{\textnormal{AUPRC}_{w}}$.
The performance is measured independently per sub-hierarchy.

Finally, Figure 8 illustrates the prediction performance measured with the
$\overline{\textnormal{AUPRC}_{w}}$, denoted as macro weighted. This score
weights the individual performance of each function according to the number of
genes associated to it. Thereby, the leaves and deeper functions in a sub-
hierarchy always get lower weight than the others. Note that the deeper a
functions is in a sub-hierarchy, the lower the predicted probabilities
becomes. The global method outperforms the locals again. The conclusion is
consistent with the three metrics, using clustering techniques to extract
features from the GCN and considering the hierarchical structure of the
biological processes seems to be key for the gene function production task.

Root | Total | Filtered
---|---|---
GO:0050896 | 239 | 124
GO:0051179 | 479 | 263
GO:0065007 | 713 | 402
GO:0008152 | 1812 | 796
GO:0009987 | 1813 | 853
Table 2: Number of extracted and filtered features used for the global method
per sub-hierarchy. Recall that the extracted features are filtered using the
mean SHAP values to select the more important with a cutoff defined by the
input constant $c$.

It has been shown in [28] that the new features built from the GCN, and the
associations between genes and functions with the spectral clustering
algorithm are key to improve the prediction performance in the gene annotation
problem (w.r.t. other features of the GCN and gene functional information).
However, the feature extraction approach presented in Section 3 produces two
different sets of features, namely, $J_{G}$ and $J_{F}$, that are combined and
used for prediction. The individual relevance of each set of features for the
gene annotation problem is analyzed by (i) looking at the distribution of the
filtered features for the global method and (ii) comparing the performance of
the prediction task using each set of features independently. Table 2 presents
the number of extracted and filtered features used for the global method per
sub-hierarchy. Recall that the features are filtered using the mean SHAP
values to select the more important ones with a cutoff defined by the input
constant $c$.

Figure 9: Distribution of the filtered features from $J_{G}$ and $J_{F}$ for
the global method per sub-hierarchy.

Figure 9 illustrates the distribution of the filtered features for the global
method per sub-hierarchy. Note that, even though the features from the
affinity graph $F$ (i.e., $J_{F}$) are more important, features from the GCN
$G$ (i.e., $J_{G}$) are also selected for all sub-hierarchies. Figure 10 shows
the prediction performance of the global HMC method trained using the features
$J_{G}$ and $J_{F}$ independently, and the proposed approach (i.e., their
combination) measured with AU($\overline{\textnormal{PRC}}$) and
$\overline{\textnormal{AUPRC}}$. The combination of both sets of features,
extracted from the GCN and the affinity graph is key to improve the
performance of the proposed approach for all sub-hierarchies.

Figure 10: Prediction performance of the global method trained using the
features $J_{G}$ and $J_{F}$ independently, and the proposed approach (i.e.,
their combination) measured with AU($\overline{\textnormal{PRC}}$) and
$\overline{\textnormal{AUPRC}}$. The performance is measured independently per
sub-hierarchy.

##  6 Related Work and Concluding Remarks

###  6.1 Related Work

Zhou et al. [41] presented an approach to predict functions of maize proteins
using graph convolutional networks. In particular, an amino acid sequence of
proteins and the GO hierarchy were used to predict functions of proteins with
a deep graph convolutional network model (DeepGOA). Their results showed that
DeepGOA is a powerful tool to integrate amino acid data and the GO structure
to accurately annotate proteins. Similarly, the work presented in [6] aims to
predict the phenotypes and functions associated to maize genes using: (i)
hierarchical clustering based on datasets of transcriptome (set of molecules
produced in transcription) and metabolome (set of metabolites found within an
organism); and (ii) GO enrichment analyses. Their results showed that
profiling individual plants is a promising experimental design for narrowing
down the lab-field gap. Gligorijević et al. [9] proposed a network fusion
method based on multimodal deep autoencoders to extract high-level features of
proteins from multiple interaction networks. This method, called deepNF,
relied on a deep learning technique that captures relevant protein features
from different complex, non-linear interaction networks. Their results showed
that extracting new features from biological networks is key to annotate gene
with functions. The work in [40] is also closely related. They presented Gene
Ontology hierarchy preserving hashing (HPHash), a gene function prediction
method that retains the hierarchical order between GO terms. It used a
hierarchy preserving hashing technique based on the taxonomic similarity
between terms to capture the GO hierarchy. Hashing functions were used to
compress the gene-term association matrix, where the semantic similarity
between genes was used to predict the functions of the genes. Their results
showed that HPHash preserves the GO hierarchy and improves prediction
performance.

In addition, the authors in [3] presented iFeature, a Python-based toolkit for
generating numerical feature representation schemes from protein sequences. It
integrated algorithms for feature clustering, selection, and dimensionality
reduction to facilitate training, analysis, and benchmarking of machine-
learning models. In a related way, Mu et al. [19] showed that feature
extraction of protein sequences is helpful for prediction of protein functions
or interactions. They introduced FEGS (Feature Extraction based on Graphical
and Statistical features), a novel feature extraction model for protein
sequences that combines graphical and statistical features. Their results
showed that similarity analysis of protein sequences has applications in the
study of gene annotation, gene function prediction, identification and
construction of gene families, and gene discovery.

###  6.2 Concluding Remarks and Future Work

By combining network-based modeling, cluster analysis, interpretable machine
learning, and hierarchical multi-label classification, the approach presented
in this paper introduces a novel method to address the gene function
prediction problem. It aims to predict the association probability between
each gene and function by taking advantage of the GCN spectral decomposition,
the information available of associations between genes and functions, and the
ancestral relations between the functions (i.e., the GO hierarchy).

A case study on Zea mayz (maize) is presented. Using the structural
information of the gene co-expression network (extracted by a spectral
clustering algorithm) and considering the hierarchical structure of the
biological processes (using HMC) seems to be the key for the improved
performance of the proposed approach. More precisely, the global HMC method,
which considers all features available for a sub-hierarchy to build a single
classifier, outperforms the other methods in relation to the three metrics
that were used (namely, AU($\overline{\textnormal{PRC}}$) ,
$\overline{\textnormal{AUPRC}}$ , and $\overline{\textnormal{AUPRC}_{w}}$).

The results presented in [28] show that the features extracted from the GCN
using spectral clustering lead to better prediction performance in the gene
function prediction task (addressed as an independent binary classification
problem per function). In this work, it has been shown that considering the
ancestral relations between functions to produce an outcome that satisfies its
hierarchical structure (i.e., complies the true-path rule or hierarchical
constraint), based on the features extracted from the GCN, improves the
performance in the gene function prediction task (addressed as a hierarchical
multi-label classification problem).

Two main lines of work can be considered for future work. First, applying the
proposed approach to identify genes associated to specific stresses (e.g., low
temperature, salinity) can help to reduce the set of candidate genes that
respond to treatments for in vivo validation. Second, exploring transfer
learning techniques (especially, domain adaptation) to enrich the building of
the classifiers using information from other organisms (datasets), not only
can lead to higher prediction performance, but also can enable the proposed
approach on organisms without a wealth of significant functional information.

## Abbreviations

AUPRC:

Area Under Precision-Recall Curve

DAG:

Directed Acyclic Graph

GO:

Gene Ontology

GCN:

Gene Co-expression Network

HMC:

Hierarchical Multi-label Classification

lcl:

Local Classifier per Level

lcn:

Local Classifier per Node

lcpn:

Local Classifier per Parent Node

## Availability of data and materials

The datasets analyzed for the current study are publicly available from
different sources. They can be found in the following locations:

  * •

Gene co-expression data of Oryza sativa Japonica is available on ATTED-II
[21].

  * •

Functional data of rice genes is available on the DAVID Bioinformatics
Resources [10].

  * •

Hierarchical data of Gene Ontology terms are available on the GOATOOLS Python
library [13].

The data collected, cleaned, and processed from the above sources as used in
the case study can be requested to the authors.

A workflow implementation is publicly available:

  * •

Project name: clustering_hmc

  * •

Project home page: https://github.com/migueleci/clustering_hmc

  * •

Operating system(s): platform independent.

  * •

Programming language: Python 3.

  * •

Other requirements: None.

  * •

License: GNU GPL v3.

## Funding

This work was partially funded by the OMICAS program: Optimización Multiescala
In-silico de Cultivos Agrícolas Sostenibles (Infraestructura y Validación en
Arroz y Caña de Azúcar), anchored at the Pontificia Universidad Javeriana in
Cali and funded within the Colombian Scientific Ecosystem by The World Bank,
the Colombian Ministry of Science, Technology and Innovation, the Colombian
Ministry of Education and the Colombian Ministry of Industry and Turism, and
ICETEX, under GRANT ID: FP44842-217-2018. The second author was partially
supported by Fundación CeiBA.

## References

  * [1] M. Ashburner, C. A. Ball, J. A. Blake, D. Botstein, H. Butler, J. M. Cherry, A. P. Davis, K. Dolinski, S. S. Dwight, J. T. Eppig, M. A. Harris, D. P. Hill, L. Issel-Tarver, A. Kasarskis, S. Lewis, J. C. Matese, J. E. Richardson, M. Ringwald, G. M. Rubin, and G. Sherlock.  Gene Ontology: Tool for the unification of biology.  Nature Genetics, 25(1):25–29, May 2000.
  * [2] S. Bhagat, G. Cormode, and S. Muthukrishnan.  Node Classification in Social Networks.  In C. C. Aggarwal, editor, Social Network Data Analytics, pages 115–148. Springer US, Boston, MA, 2011.
  * [3] Z. Chen, P. Zhao, F. Li, A. Leier, T. T. Marquez-Lago, Y. Wang, G. I. Webb, A. I. Smith, R. J. Daly, K.-C. Chou, and J. Song.  iFeature: A Python package and web server for features extraction and selection from protein and peptide sequences.  Bioinformatics, 34(14):2499–2502, July 2018.
  * [4] H. Cho, B. Berger, and J. Peng.  Diffusion Component Analysis: Unraveling Functional Topology in Biological Networks.  In RECOMB 2015, pages 62–64. Springer, Cham, Apr. 2015.
  * [5] H. Cho, B. Berger, and J. Peng.  Compact Integration of Multi-Network Topology for Functional Analysis of Genes.  Cell Systems, 3(6):540–548.e5, Dec. 2016.
  * [6] D. F. Cruz, S. De Meyer, J. Ampe, H. Sprenger, D. Herman, T. Van Hautegem, J. De Block, D. Inzé, H. Nelissen, and S. Maere.  Using single-plant-omics in the field to link maize genes to functions and phenotypes.  Molecular Systems Biology, 16(12), Dec. 2020.
  * [7] M. Deng, K. Zhang, S. Mehta, T. Chen, and F. Sun.  Prediction of Protein Function Using Protein–Protein Interaction Data.  Journal of Computational Biology, 10(6):947–960, Dec. 2003.
  * [8] Gene Ontology Consortium.  The Gene Ontology Resource: 20 years and still GOing strong.  Nucleic Acids Research, 47(D1):D330–D338, Jan. 2019.
  * [9] V. Gligorijević, M. Barot, and R. Bonneau.  deepNF: Deep network fusion for protein function prediction.  Bioinformatics, 34(22):3873–3881, 2018.
  * [10] D. W. Huang, B. T. Sherman, and R. A. Lempicki.  Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources.  Nature Protocols, 4(1):44–57, Jan. 2009.
  * [11] H. Jia, S. Ding, X. Xu, and R. Nie.  The latest research progress on spectral clustering.  Neural Computing and Applications, 24(7-8):1477–1486, June 2014\.
  * [12] S. S. Khan and M. G. Madden.  A Survey of Recent Trends in One Class Classification.  In L. Coyle and J. Freyne, editors, Artificial Intelligence and Cognitive Science, volume 6206, pages 188–197. Springer Berlin Heidelberg, Berlin, Heidelberg, 2010.
  * [13] D. V. Klopfenstein, L. Zhang, B. S. Pedersen, F. Ramírez, A. Warwick Vesztrocy, A. Naldi, C. J. Mungall, J. M. Yunes, O. Botvinnik, M. Weigel, W. Dampier, C. Dessimoz, P. Flick, and H. Tang.  GOATOOLS: A Python library for Gene Ontology analyses.  Scientific Reports, 8(1):10872, Dec. 2018.
  * [14] J. Levatić, D. Kocev, and S. Džeroski.  The importance of the label hierarchy in hierarchical multi-label classification.  Journal of Intelligent Information Systems, 45(2):247–271, Oct. 2015.
  * [15] S. Lundberg and S.-I. Lee.  A Unified Approach to Interpreting Model Predictions.  arXiv:1705.07874 [cs, stat], Nov. 2017.
  * [16] S. M. Lundberg, G. Erion, H. Chen, A. DeGrave, J. M. Prutkin, B. Nair, R. Katz, J. Himmelfarb, N. Bansal, and S.-I. Lee.  From local explanations to global understanding with explainable ai for trees.  Nature Machine Intelligence, 2(1):2522–5839, 2020.
  * [17] F. Luo, Y. Yang, J. Zhong, H. Gao, L. Khan, D. K. Thompson, and J. Zhou.  Constructing gene co-expression networks and predicting functions of unknown genes by random matrix theory.  BMC Bioinformatics, 8(1):299, Dec. 2007.
  * [18] P. Mills.  Solving for multi-class: A survey and synthesis.  arXiv:1809.05929 [cs, stat], Jan. 2021.
  * [19] Z. Mu, T. Yu, X. Liu, H. Zheng, L. Wei, and J. Liu.  FEGS: A novel feature extraction model for protein sequences and its applications.  BMC Bioinformatics, 22(1):297, Dec. 2021.
  * [20] N. Murugesan, I. Cho, and C. Tortora.  Benchmarking in Cluster Analysis: A Study on Spectral Clustering, DBSCAN, and K-Means.  In Data Analysis and Rationality in a Complex World, pages 175–185. Springer, Cham, 2021.
  * [21] T. Obayashi, Y. Aoki, S. Tadaka, Y. Kagaya, and K. Kinoshita.  ATTED-II in 2018: A Plant Coexpression Database Based on Investigation of the Statistical Property of the Mutual Rank Index.  Plant and Cell Physiology, 59(1):e3–e3, Jan. 2018.
  * [22] T. Obayashi and K. Kinoshita.  COXPRESdb: A database to compare gene coexpression in seven model animals.  Nucleic Acids Research, 39(Database):D1016–D1022, Jan. 2011.
  * [23] M. Oti, J. van Reeuwijk, M. A. Huynen, and H. G. Brunner.  Conserved co-expression for candidate disease gene prioritization.  BMC Bioinformatics, 9(1):208, 2008.
  * [24] G. A. Petsko.  Guilt by association.  Genome Biology, 10(4):104, 2009.
  * [25] M. Z. Rodriguez, C. H. Comin, D. Casanova, O. M. Bruno, D. R. Amancio, L. d. F. Costa, and F. A. Rodrigues.  Clustering algorithms: A comparative approach.  PLOS ONE, 14(1):e0210236, Jan. 2019.
  * [26] M. Romero, J. Finke, M. Quimbaya, and C. Rocha.  In-silico Gene Annotation Prediction Using the Co-expression Network Structure.  In Complex Networks and Their Applications VIII, pages 802–812. Springer, Cham, 2020.
  * [27] M. Romero, J. Finke, and C. Rocha.  A top-down supervised learning approach to hierarchical multi-label classification in networks.  Applied Network Science, 7(1):8, Dec. 2022.
  * [28] M. Romero, Ó. Ramírez, J. Finke, and C. Rocha.  Supervised Gene Function Prediction Using Spectral Clustering on Gene Co-expression Networks.  In R. M. Benito, C. Cherifi, H. Cherifi, E. Moro, L. M. Rocha, and M. Sales-Pardo, editors, Complex Networks & Their Applications X, volume 1016, pages 652–663. Springer International Publishing, Cham, 2022\.
  * [29] A. G. Rust, E. Mongin, and E. Birney.  Genome Annotation Techniques: New Approaches and Challenges.  Drug Discovery Today, 7(11):S70–S76, May 2002.
  * [30] C. N. Silla and A. A. Freitas.  A survey of hierarchical classification across different application domains.  Data Mining and Knowledge Discovery, 22(1-2):31–72, Jan. 2011.
  * [31] J. M. Stuart.  A Gene-Coexpression Network for Global Discovery of Conserved Genetic Modules.  Science, 302(5643):249–255, Oct. 2003.
  * [32] G. Valentini.  True Path Rule Hierarchical Ensembles.  In Multiple Classifier Systems, pages 232–241. Springer, Berlin, Heidelberg, 2009.
  * [33] S. van Dam, U. Võsa, A. van der Graaf, L. Franke, and J. P. de Magalhães.  Gene co-expression analysis for functional classification and gene–disease predictions.  Briefings in Bioinformatics, page bbw139, Jan. 2017.
  * [34] K. Vandepoele, M. Quimbaya, T. Casneuf, L. De Veylder, and Y. Van de Peer.  Unraveling Transcriptional Control in Arabidopsis Using cis-Regulatory Elements and Coexpression Networks.  Plant Physiology, 150(2):535–546, June 2009.
  * [35] C. Vens, J. Struyf, L. Schietgat, S. Džeroski, and H. Blockeel.  Decision trees for hierarchical multi-label classification.  Machine Learning, 73(2):185–214, Nov. 2008.
  * [36] D. Xu, Y. Shi, I. W. Tsang, Y.-S. Ong, C. Gong, and X. Shen.  Survey on multi-output learning.  IEEE Transactions on Neural Networks and Learning Systems, 31(7):2409–2429, 2020.
  * [37] M. Yandell and D. Ence.  A beginner’s guide to eukaryotic genome annotation.  Nature Reviews Genetics, 13(5):329–342, May 2012.
  * [38] S. Yon Rhee, V. Wood, K. Dolinski, and S. Draghici.  Use and misuse of the gene ontology annotations.  Nature Reviews Genetics, 9(7):509–515, July 2008.
  * [39] Yu and Shi.  Multiclass spectral clustering.  In Proceedings Ninth IEEE International Conference on Computer Vision, pages 313–319 vol.1, Nice, France, 2003. IEEE.
  * [40] Y. Zhao, G. Fu, J. Wang, M. Guo, and G. Yu.  Gene function prediction based on Gene Ontology Hierarchy Preserving Hashing.  Genomics, 111(3):334–342, 2019.
  * [41] G. Zhou, J. Wang, X. Zhang, M. Guo, and G. Yu.  Predicting functions of maize proteins using graph convolutional network.  BMC Bioinformatics, 21(S16):420, Dec. 2020.
  * [42] Y. Zhou, J. A. Young, A. Santrosyan, K. Chen, S. F. Yan, and E. A. Winzeler.  In silico gene function prediction using ontology-based pattern identification.  Bioinformatics, 21(7):1237–1245, Apr. 2005.