# PhenoLinker: Phenotype-Gene Link Prediction and Explanation using
Heterogeneous Graph Neural Networks

(  Jose L. Mellina Andreu1, Luis Bernal1, Antonio F. Skarmeta1, Mina Ryten2,3,
Sara Álvarez4, Alejandro Cisterna García1*, Juan A. Botía1,3*† 1Departamento
de Ingeniería de la Información y las Comunicaciones, Universidad de Murcia,
Murcia, Spain. 2Genetics and Genomic Medicine, Great Ormond Street Institute
of Child Health, University College London, London, WC1E 6BT, UK. 3Department
of Neurodegenerative Diseases, UCL Institute of Neurology, London WC1E 6BT, UK
4NIMGenetics, Genómica y Medicina, SL, Madrid, Spain * The authors contributed
equally †Corresponding author at: Departamento de Ingeniería de la Información
y las Comunicaciones, Universidad de Murcia, Murcia, Spain. E-mail address:
<EMAIL_ADDRESS>(J.A. Botía). )

###### Abstract

The association of a given human phenotype to a genetic variant remains a
critical challenge for biology. We present a novel system called PhenoLinker
capable of associating a score to a phenotype-gene relationship by using
heterogeneous information networks and a convolutional neural network-based
model for graphs, which can provide an explanation for the predictions. This
system can aid in the discovery of new associations and in the understanding
of the consequences of human genetic variation.

##  1 Introduction

In the era of artificial intelligence (AI), deep learning has led to
remarkable advancements in various fields, with examples such as AlphaFold
[1], GPT-4 [2] or the very recent AlphaMissense [3]. In all these advances,
all input data had a tabular nature, i.e., all training examples were
expressed as a regular, component-based, vector of values. However, as we
confront modeling problems where some of the central data structures go beyond
that static and rigid structure to become graphs, the need for specialized
techniques emerges. Graph Neural Networks (GNNs) have been advancing rapidly
in recent years, improving the results of classical deep learning techniques
by combining the use of tabular data with graph structures, processing data
that can be represented as graphs, such as social or citation networks. In the
field of biology, these models are rapidly improving results in a multitude of
tasks, such as the prediction of unlabeled proteins or the generation of
synthetic molecules through actual molecular graph learning [4], prediction of
protein-phenotype associations with protein-protein associations networks [5],
or even for drug-target interactions [6].

We are concerned with the task of predicting new associations between
phenotypes and genes. Phenotypes are observable traits above the molecular
level, such as anatomy and behavior [7]. In our context, a gene-phenotype
association represents a critical link in understanding genetic diseases, as
it denotes that pathogenic mutations within the gene in question hold the
potential to give rise to the associated phenotype [8, 9]. The impact of
genetics on human health is broad and challenging to measure primarily because
it encompasses various types of diseases caused either entirely or partially
by genetic factors such as complex disorders, chromosomal disorders or single-
gene disorders [10, 11]. For instance, as indicated by certain studies that
concentrate solely on rare diseases, approximately 71.9% of which have a
genetic basis, the population prevalence of rare diseases is estimated to be
around 3.5–5.9% [12]. Hence, focusing exclusively on genetic rare diseases, we
can estimate that the population prevalence would be in the range of
approximately 2.5-4.2%. In the realm of genetic diagnostics using whole-exome
sequencing information, identifying gene-phenotype associations becomes a
crucial step [13, 14, 15]. Here, clinicians, prompted by a patient’s symptoms
and phenotypes, curate a list of potential phenotypes, which, through
previously documented gene-phenotype associations, ultimately lead to
candidate genes. These genes then lead to causative variants or mutations.
Therefore, phenotypes lead to genes and genes lead to possible candidate
mutations from a comprehensive and carefully manually curated catalog. And a
successful diagnostic happens if the geneticist manages to find a causal
mutation given that the patient has that mutation and manifest the phenotypes
associated with it, through the given genes. Therefore, the importance of
comprehensive gene-phenotype associations becomes evident in the diagnostic
process. When our knowledge about what we know of gene-phenotype associations
in regard to the target disease is incomplete, due to unavailability or lack
of reliable documentation, we fail to provide the requisite diagnostic
information. As a precise diagnosis leads to a more specific and maybe
personalized treatment, a successful diagnosis is paramount. Unknown
phenotype-gene associations contribute as one of the factors responsible for
leaving undiagnosed approximately 60% of the patients under suspicion of
suffering a genetic disease [16]. On top of that, it is also essential, when
looking for the causative mutation, to keep the number of genes to inspect
during the process (i.e., the so-called gene panel for the phenotypes
implicated), in a specific set. Expanding the set of genes with the aim of not
excluding any possible causal mutation will result in a broader range of
potential mutations. While increasing that probability marginally, the
increase of work overload for the expert geneticist may be unacceptable.
Consequently, achieving a balance between precision and recall proves pivotal
in this undertaking [17]. Therefore, there is a need for more complete yet
highly reliable gene-phenotype associations to increase the yield of genetic
diagnosis without injecting an additional burden to medical experts.

To formally model what is relevant to describe associations between phenotypes
and genes, we have to consider the Human Phenotype Ontology (HPO) [18]. HPO is
an ontology whose structure is a directed acyclic graph. It was created to
describe the spectrum of all possible human phenotypic abnormalities. It is
based on a standardized vocabulary to provide a unified representation of
phenotypes and their relationships, mainly the is-a relation, that is used to
express that phenotype a is a subtype of phenotype b if we say b is-a a. We
have also considered as potentially useful other ontology-based sources of
information about molecular processes and diseases, such as the Gene Ontology
(GO) [19, 20], the Kyoto Encyclopedia of Genes and Genomes (KEGG) [21] or the
Human Disease Ontology (DO) [22]. However, none of these sources include gene-
phenotype associations suitable for genetic diagnosis that we can use to
predict new gene-phenotype associations. Based on HPO, several gene-phenotype
prediction methods have been proposed in the last few years, such as HPODNets
[5], CADA [23] and others [24, 25, 26].

Here, we have designed a totally new approach to predict gene-phenotype
associations by considering Heterogeneous Information Networks (HINs), so that
we have a network with nodes of different classes and with different
attributes, such as phenotypes and genes, based on information from HPO. HINs
are extremely interesting in this case for a number of reasons: i) Having
phenotypes and genes in the same network allows treating the gene-phenotype
associations seamlessly within the GNN. ii) The possibility of describing each
node in the network with specific attributes leads the door open to a more
detailed description of genes and phenotypes and hence, to a more
comprehensible set of gene-phenotype associations and iii) the GNN training
algorithm is capable of capturing how specific gene features impact the
potential gene-phenotype associations. This new layer of information is
beneficial for obtaining more accurate predictions, and can contribute to
reach an understanding, through explainable AI, of which features drive a gene
to be associated with a phenotype. Following these hypotheses, we developed a
model called PhenoLinker based on heterogeneous structure graph neural
networks (HSGNNs) to make the most out of this HIN (see Fig. 1). Given a pair
gene-phenotype used as input, the HSGNN developed in this research provides a
quality score and explanations for that pair that sheds light on how the
prediction was made by the neural network. The predictions of the model,
scores, explanations and different visualization approaches on a certain
dataset are made available to the research community in a web application we
have developed via Hugging Face Spaces.

Figure 1: Overview of the overall approach for PhenoLinker. A graph is created
from the HPO ontology structure which describes all phenotypes in the form of
a hierarchy and includes known genes associated to phenotypes. Starting from
this structure, we decorate the corresponding gene nodes with extended
features gathered from different sources reflecting different aspects of their
function: genomics, gene expression and genetics. From the resulting graph, a
HIN resembling that graph is created. That HIN is exploited by the HSGNN of
step 2 to encode the subgraph of interest for the phenotype and the gene into
vectors of real numbers (i.e., embeddings). This HSGNN was optimized to result
in a SAGE based (see methods) three layer codecs for phenotype and gene. With
the embeddings for both phenotype and gene, the HSGNN is then capable of
predicting whether the input pair phenotype-gene should actually exist or not:
the score for a given link is calculated as the sigmoid function of the dot
product of gene and phenotype embeddings. Finally, the prediction can be
explained on the basis of which attributes and their values were used to
describe the phenotype and the gene, with a score for each attribute
reflecting its importance in the whole process.

##  2 Materials

We use the December 15, 2022 version of the HPO ontology, directly available
on their website, with the ontology file (hp.obo) and the file containing the
gene relationships (phenotype_to_genes.txt). For gene attributes, information
is gathered from several sources as it is done in [27]. We have gene-specific
variant frequency information belonging to the Gnomad database [28],
information of gene and transcript structural complexity from GENCODE (Ensembl
version 72) [29] and the HEXEvent database [30], gene-specific protein-protein
interactions information from the STRING database [31] and tissue-specific
gene expression and co-expression information from the GTEx V6 gene expression
dataset. For a detailed view of the attributes see Supplementary Table S1.

A BioBERT model [32] available from Hugging Face is used to create the
phenotype attributes. This model was previously trained with the SNLI, MNLI,
SCINLI, SCITAIL, MEDNLI, and STSB datasets, and creates 768-dimensional
embeddings that we use as phenotype attributes from their textual description.
The interactive online application is based on the Genomics England amber
panel Intellectual disability - microarray and sequencing v5.204 [33].

For the validation of our model, we have also accessed data from the
Gene2Phenotype [34] database, using the relationships of all the panels that
were not included in the used HPO version. This data was downloaded on
September 20, 2023, as the data is generated on the fly with the most recent
information directly from the website.

##  3 Methods

In this section, we describe the components of our model and how we applied
it. We introduce the main architecture of our model, explaining the process of
creating the dataset and the neural network architecture that we use to make
the predictions, as well as the additional layer that implements the
interpretation of the results.

The fundamental aim that drives our research is the discovery of new
undiscovered associations between genes and phenotypes. To enable this using
supervised machine learning (ML), we turn the original problem into a link
prediction task on a heterogeneous graph, so that we are interested in links
not belonging to the graph that our model predicts as belonging to it (i.e.,
the ML model predicts a high quality for that link, and subsequently, for the
corresponding pair gene-phenotype). To construct the model for the link
prediction task, the necessary steps are as follows: (1) creation of the graph
with phenotypes and genes emerging from HPO, the HIN, (2) creation of the
HSGNN model for heterogeneous graphs capable of generating suitable embeddings
for genes and phenotypes, (3) inference of scores for all possible and
potentially new pairs gene-phenotype and (4) generation of interpretation of
the most promising predictions.

###  3.1 HIN graph of phenotypes and genes generation

We use HPO as the basis for our graph construction as it can be seen in the
upper left part of Fig. 1 at Step 1. The HPO ontology stores information about
phenotypes as an ontology, which can be viewed as an implicit knowledge graph
with the is-a relation. We extract this knowledge graph and use it as the base
graph to which we add gene nodes and edges between genes and phenotypes using
the gene annotations provided also by HPO, where they incorporate in tabular
form the existing phenotype-gene relationships.

Apart from the is-a relation between the individuals (in this case, the
phenotypes), the HPO ontology also contains attributes describing the objects.
Specifically, the ontology stores a textual description for each phenotype as
an attribute, which we retrieve and use as the attribute of the phenotype
nodes after transforming them into an embedding. For example, for the ontology
term HP:0001945 (Fever), the ontology includes the description Body
temperature elevated above the normal range.

To transform those descriptions into semantics preserving vectors of real
values, we used a pre-trained BioBERT [32] model. This natural language
processing model transforms English language sentences of variable length into
a 768-dimensional embedding. Since this BioBERT is pre-trained on biomedical
texts, we use it without any further fine-tuning. We use the embedding as the
vector of attributes of the node for the phenotype. This allows the
incorporation of numerical attributes without having to generate a specific
model or fine-tuning it.

With respect to gene nodes at the HIN, and which attributes we use to describe
them we have compiled the available information and generated 64 attributes
divided into 5 groups: Gnomad parameters on the genetic burden conveyed by the
gene (4), gene expression across a variety of human control tissues (47
binary), parameters about the complexity of the genome locus that harbors the
gene (2), the average level of interaction of the gene when translated into
proteins (1) and Genomic Information (10). We designed this attribute set with
the goal of indirectly describing gene function in relation with the
phenotypes implied. All variables were scaled and centered before being
introduced into the graph.

##  4 Embedding generation for genes and phenotypes

After we generate the graph that holds all the known gene-phenotype
associations, with their corresponding attribute values, we need to
appropriately learn how to encode the latent structures based on the graph
that are used to encode genes and phenotypes into vectors of real numbers. In
this approach, producing a score for any pair gene-phenotype is based on the
cosine similarity of the corresponding embedding-based codification of both
gene and phenotype. Therefore, the next step is to design codecs for genes and
phenotypes that capture both the information provided through the graph
structure and through the attributes of all nodes considered in their
codification.

For this, we propose an HSGNN architecture based on SAGEConv [35] layers,
which is implemented in Pytorch Geometric using a new data structure called
HeteroData, that allows us to create different sets of nodes with different
attributes.

The SAGEConv operator uses a set of aggregator functions that learn to
aggregate feature information from a node’s local neighborhood in the graph.
Given a certain node $v$ in a graph $G=(V,E)$ we define the neighborhood as
$N(v)=\\{u\in V|(v,u)\in E\\}$, where $V$ and $E$ are the sets of nodes and
edges in the graph, respectively. The node $v$ has a certain vector of numeric
attributes $x_{v}$ that represents the attributes of the node. The operation
of the SAGEConv operator is defined as

$x_{v}^{l+1}\leftarrow W_{1}^{l}\cdot
x_{v}^{l}+W_{2}^{l}\cdot\text{mean}_{j\in\mathcal{N}(v)}x_{j}^{l}+\text{bias}^{l},$

where $l$ is the layer, $W1,W2$ are weight matrices and $\text{bias}^{l}$ is
the bias of the linear operation. The dimensions of each matrix and the bias
are determined by the input and output dimensions of the particular layer.
Thus, for an input of size $n$ and an output of size $m$, the matrices $W1,W2$
belong to $M_{n\times m}$, and the bias is a vector of $m$ dimensions.

In our architecture, we define the neighborhood in a heterogeneous graph
depending on the specific subgroup of edges. Specifically, for a phenotype
node, we have two subgroups of edges: those that connect the node with another
phenotype and those that connect it with a gene. To handle these subgroups, we
have a separate SAGEConv operator for each kind of edge. The final embedding
of the phenotype nodes is generated by adding the embeddings from the layers.
This means that the resulting embedding for the phenotype encodes both
information from the phenotype itself and from the associated genes.

After each layer of SAGEConv operators, a LeakyReLU [36] activation function
is selected. The final embedding is directly generated by a final
heterogeneous SAGEConv layer. The final structure can be seen in Fig. 1. By
empirical experimentation, it has been decided to use 3 intermediate SAGE
layers followed by a final one that generates the embeddings, and 64
dimensions for the hidden channels of the model. This is shown in the results
section.

###  4.1 Prediction generation

The final step involves predicting a score for a given edge in the
heterogeneous graph. To achieve this, we leverage the node embeddings
generated in the previous step and calculate the dot product between the
embeddings. Subsequently, we use a sigmoid function to produce a probability
for the given pair.

Formally, given a phenotype p with a representation $(p_{1},\dots,p_{n})$ and
a gene $g$ with a representation $(g_{1},\dots,g_{n})$, the model computes the
dot product and then is transformed into a probability using the sigmoid
function as

$\sigma(p\cdot g)=\frac{1}{1+e^{-\sum_{i=1}^{n}}p_{i}g_{i}}.$

###  4.2 Training and evaluating the model

The task can be viewed as a binary classification problem for a given edge in
the network. The positive class corresponds to the edges that belong to the
network, and the negative class to the edges not present in the network. We
use as the loss function the binary cross-entropy for each edge in the
training set. Formally, for a certain vector of predictions
$\hat{Y}=(\hat{y_{1}},\dots,\hat{y_{n}})$ and the corresponding reference
values Y=(y1, … , yn), the objective function of optimization is

$\frac{1}{n}\sum_{i=1}^{n}(y_{i}\log(\hat{y_{i}})+(1-y_{i})\log(1-\hat{y_{i}})).$

This problem has some particularities derived from the fact that the potential
set of negative examples (i.e., the whole set of links not present in the HPO-
derived graph of phenotypes and genes), indeed has unknown positive examples
that we want to detect. In order to create a training dataset, we need
positive cases (edges that are in the network) and negative cases (edges that
are not). Positive cases are randomly selected from the existing links. But
for the negative cases the situation is slightly more complex as in this case
we focus on phenotype-gene pairs which are not connected and, from those, the
absent links are used.

As usual in supervised experiments, we need to organize the training data,
i.e. the edges, into training, validation and test sets. This split is done
using RandomLinkSplit from Pytorch Geometric, where an important parameter is
the negative sampling ratio to include in the sets. This is the percentage of
the total of potentially usable links not existing in the graph, that we may
use as negative examples. Each training example in the training data is
composed of a series of message passing edges and one supervision edge.
Message passing edges are used by the GNN to pass the message. They are the
edges through which the network passes the information in each layer. And
supervision edges are used in the loss function for backpropagation. During
training, message-passing edges are used to transmit the information. This
means that these are the edges that will be incorporated into the SAGE layers
to generate the embeddings in the model. Once the embeddings are generated,
the supervision edges are used to calculate the training error and therefore,
they are used to calculate the error updates that propagate backwards through
the HSGNN weights. For validation, we use all the previous edges as input to
the model to calculate the embeddings, and we have new validation edges that
are the ones that will give us the metrics.

##  5 Explanation of link predictions

After we use the HSGNN for scoring a gene-phenotype link of interest, the
model is capable of explaining the produced score in terms of the particular
contribution of each node attribute, both for phenotypes and genes. This
useful feature we have incorporated in our model is based on the Integrated
Gradients [37] algorithm, available from the Captum [38] library through
Pytorch Geometric. This algorithm returns a weight associated with each
attribute in reference to its importance when scoring the link between the
gene and the phenotype.

For the sake of usability of our predictions, these interpretations are
displayed in an interactive online application using Hugging Face Spaces,
where a 3D representation of the embeddings is shown using the t-SNE [39]
method, and in addition, the 10 most important attributes are displayed with
their scores for a number of genes and phenotypes. It can be accessed at
https://huggingface.co/spaces/Melland/PhenoLinker.

##  6 Results

###  6.1 Data preparation

We used Pytorch Geometric 2.4.0 as the basis for our architecture, as it is a
popular library for graph neural networks in Python. To train the neural
network, we used a cross-entropy loss function, Adam optimizer, a learning
ratio of 0.001, 15 epochs and a batch size of 2048. The core of our training
data is the phenotypes within the HPO ontology in its December 15, 2022
version and the genes of the HPO annotations of the same version. Only the
genes with no missing values in the attributes are used. Our final dataset has
16810 phenotypes, 4619 genes, 42142 relationships between phenotypes and
721540 relationships between genes and phenotypes. We define negative sampling
ratio as the ratio of examples of the negative class to the positive class.
Given that, we have chosen a negative sampling ratio of 4 when building the
negative samples, implying 20% of positive class examples in our training set
for the link prediction task, and 30% of supervision edges in the training
set. Our estimate is that in the full model there is a negative ratio of 107.
Therefore, the learning process is clearly imbalanced. To alleviate this, we
subsample the negative class so the negative ratio is only 4. This choice
responds to a trade-off between balancing the classes and somewhat respecting
the imbalanced nature of the problem.

To estimate the predictive power of our model, we used the area under the
precision-recall curve (AUCPR). This curve has been widely used as a preferred
metric for evaluating link prediction models in large graphs [40]. The link
prediction task suffers from the problem of the existence of mislabeled edges
(false negatives and false positives in the classification task), since a
relationship could be undiscovered but should exist in the graph, or the
relationship could be wrongly added. The use of the AUCPR metric partially
alleviates the problem for the false negatives and is therefore preferred in
these tasks. The Average Precision (AP) formula is used to calculate this
area. It summarizes the curve as the weighted mean of precisions achieved at
each threshold, with the increase of recall from the previous threshold used
as the weight.

We designed a hyperparameter grid of possible values with the number of SAGE
layers that generates the embeddings of the nodes between 1 and 8, as well as
with the possible dimensions of the embeddings and hidden channels of the
network between 16 and 256 (see Fig. 2, panel a). The training process showed
that the model did not suffer from overfitting to the training data (Fig. 2,
panel b). The setting was repeated 10 times with random sampling, and the mean
and standard deviation are shown. Following the results, 3 layers and 64
dimensions for the embeddings were chosen, so that we have selected the
simplest model without losing precision.

###  6.2 Prediction performance goes beyond the state of the art.

Figure 2: a) Comparative plot of model hyperparameters. The mean AUCPR of 10
random repetitions is calculated and the mean and standard deviation are taken
for each combination. We have chosen 3 layers and an embedding size of 64
dimensions. b) Plot of the evolution of the AUCPR in the training and
validation sets during 15 training epochs. The evolution of the loss function
is also shown. Training set of 80% and validation set of 20% using the
hyperparameters chosen from a)

To compare the results of our model, we used the results obtained by HPODNets
[5] in their paper. To the best of our knowledge, HPODNets is the deep
learning-based approach that is used in a similar task to ours, and seems to
be the most powerful. In the corresponding paper, its performance is compared
with that of other prediction methods. This comparison is restricted to only
the Phenotypic Abnormality (PA) branch terms of HPO whose phenotypes were
additionally filtered leaving out those with less than 10 associated proteins
through UniProt (proteins that were annotated in Swiss-Prot). With this set,
the model was trained with the complete dataset, with all possible negative
class edges (which implies a negative ratio of 26, that is six times greater
than the one our model is trained on, i.e. 4).

Two metrics were chosen to compare our results with HPODNets: AUCPR and F1
score. We use AUCPR because it is the metric we employ thoroughly in this
paper. We also use F1 because it is one of the other metrics they use. Each
was reported in two ways: (i) Macro-averaged (M): calculate metrics for each
term (for each gene), and calculate the mean; (ii) Micro-averaged (m):
vectorize the predicted score for each pair and calculate metrics based on the
resulting vector. The dataset used has been emulated as closely as possible.
We have used the December 15, 2022 version of HPO and the current UniProt and
Swiss-Prot proteins to filter out the phenotypes, so we have obtained a
dataframe with 4574 phenotypes and 4619 genes. On this dataset we have used
the model previously trained on the global dataset, with 3 layers and
64-dimensional embeddings, using all the attributes for the nodes of the
graph.

The comparison is found in Table 1, where it can be seen that our model
obtains the best results by a wide margin in M-AUCPR with 50.37%, m-AUCPR with
53.85% and M-F1 with 64.33%, but do not obtain the best result in m-F1 with
33.08% whilst HPODNets obtains a 36.15%. Since this data is really unbalanced,
with a much higher proportion of negative edges, the m-F1 of our model gets
much worse, as it has been designed to detect missing edges in the network and
it detects many of these edges as positive. The F1 metric can be misleading
for these unbalanced sets as it does not account for true negatives and false
positives are not necessarily errors in this context.

Method | M-AUCPR | m-AUCPR | M-F1 | m-F1
---|---|---|---|---
deepNF | 0.28833 | 0.24189 | 0.36646 | 0.30281
DeepMNE | 0.31354 | 0.26899 | 0.39181 | 0.32425
BIONIC | 0.30441 | 0.26587 | 0.38448 | 0.31973
LP | 0.29329 | 0.25042 | 0.37203 | 0.30403
RANKS | 0.27792 | 0.23105 | 0.36329 | 0.29188
Mashup | 0.33066 | 0.29211 | 0.40789 | 0.34948
GeneMANIA | 0.33314 | 0.23328 | 0.40995 | 0.29647
HPODNets | 0.35754 | 0.31308 | 0.43397 | 0.36152
PhenoLinker | 0.50379 | 0.53855 | 0.64326 | 0.33080
Table 1: Performance comparison. The results of HSGNN are summarized over 5
trials. Results of previous methods are extracted from [5] using a weighted
mean. The boldface items in the table represent the best performance.

###  6.3 Temporal validation performance is above state of the art approaches.

We wanted to assess whether our model is able to predict absent links that
should eventually be discovered in the future. Therefore we performed what we
call a temporal validation. For this, we adopted a strategy where the training
set is a version of HPO (December 15, 2022), and as links discovered a
posteriori, we define as test sets the new and deleted links that we find in
posterior versions of HPO. In this way, we can produce rough estimates of
whether the model is able to predict these new additions and false relations
(detected and removed by the HPO maintainers).

There are 7 versions of HPO posterior to the one we have used to train the
model, although we discarded one of them (April 5, 2023) because it presents
abnormally large differences with respect to the other versions (it discards
10% of the total relationships), and it seems to be wrongly released as this
change was fixed again in the following versions. In all the other cases, only
the phenotypes and genes belonging to the training set are considered in the
validation set. Accordingly, the links used to evaluate the model with, are
those that incorporate both the new relations (positive class) and those
eliminated (negative class) for those phenotypes and genes only.

The results obtained are shown in Table 2. It can be seen that the version
that gives the best results is the 2023-06-17 version. As time passes and the
versions are later, the differences between the versions increase, and as the
differences increase so does the performance of the model, because our model
is trained on older data. Finally, the model behaves worse with the deleted
data since the model is training with them as positively labeled edges, so
this could be considered as erroneous information. We can observe an almost
perfect positive correlation between the values of the AUCPR column and those
for the value of the ratio of added/removed links (Pearson correlation 0.9839,
P= 0.0003854).

We can compare the results with HPODNets because they also did a temporal
validation too. They used as the base for training data the release of
February 12, 2019. Their predictions were tested with the data release on
October 12, 2020, using the m-AUCPR and M-AUCPR metrics. These metrics for the
HPODNets method using the weighted mean according to term groups (they are
grouped by the number of associated proteins) is 11.5% for M-AUCPR and 8.3%
for m-AUCPR. Thus, PhenoLinker shows better results than those offered by
HPODNets (see Table 2).

YY-mm-dd | AUCPR | recall | precision | F1 | New relations | Removed relations | Ratio added/removed
---|---|---|---|---|---|---|---
2023-01-27 | 0.5906 | 0.6987 | 0.6231 | 0.6587 | 2310 | 1039 | 2,22
2023-06-06 | 0.7947 | 0.3872 | 0.7731 | 0.5160 | 42656 | 5409 | 7,88
2023-06-17 | 0.8027 | 0.3871 | 0.7806 | 0.5175 | 43059 | 5179 | 8,31
2023-07-21 | 0.7611 | 0.4365 | 0.7504 | 0.5520 | 62984 | 9880 | 6,37
2023-09-01 | 0.7449 | 0.4254 | 0.7266 | 0.5366 | 58630 | 10136 | 5,78
2023-10-09 | 0.7049 | 0.4289 | 0.6872 | 0.5281 | 60631 | 12705 | 4,77

Table 2: Temporal validation results. AUCPR, recall, precision and F1 are
shown in the sets formed with the new HPO versions. It is shown how many
relationships are added and removed with respect to the training version
(December 15, 2022), and the ratio between them.

###  6.4 Using phenotype and gene-level attributes is effective

One of the innovations introduced by this paper is the use of attributes in
the nodes of the GNN to provide a more informative description of both
phenotypes and genes. To assess whether they contribute to creating better
models, we faced the following questions: Is it any better to use attributes
at all? Are attribute values really meaningful? What is the particular
contribution of each attribute as measured separately? To try to shed some
light on the answers to these questions, respectively, we estimated the
precision and recall of a model with no attributes, a model with attributes
whose values were assigned randomly (i.e., to simulate totally uninformative
attributes) and models with a single attribute so we test how productive is
each attribute alone. We followed this approach separately for the attributes
at phenotypes (see Fig. 3, panel a) and genes (Fig. 3, panel b).

Therefore, with respect to attributes at phenotypes, to assess whether the
inclusion of the phenotype attributes coded into vectors of R768 through a
BioBERT model is really effective in PhenoLinker, we compared it with a model
with no attributes and another one with attributes with randomly assigned
values, as we can see in Fig. 3., panel a). We see that the AUCPR is
significantly higher when we include the attributes and their values are not
random. We also see that it is better not to use attributes than to use noisy
ones, as we can observe with the loss of precision using random attributes.

Figure 3: Comparison of models with different groups of attributes. The
experiment was developed from 10 randomized repetitions with a training set of
80% and test set of 20%. a) We compare the precision-recall plot of the random
model, a model with randomized attributes and the model where we add the
embeddings of the phenotypes as attributes. AP is the Average Precision
formula for the area under the Precision-Recall curve. c) and d) Graphs
showing precision and recall for the different groups of attributes. The mean
over the 10 replicates is plotted, and the standard deviation is shown as an
error bar. The distributions were compared using Welch’s t-test where we
compared the means of the distributions, and symbols representing the p-value
of significance are shown.

With respect to attributes at genes, we again compare the predictions for each
group of gene attributes with empty and random models in Fig. 3, panel b). All
of the groups of attributes show better results than the model with no
attributes, which again is better than one using random values in the
attributes. The AUCPR for the random model is 89%, for the empty model is 90%
and all models with gene attributes get an AUCPR in the interval [93-95]%.
Interestingly, the differences with respect to recall are not significant for
most groups (see Fig. 3, panel c), except for Locus, Protein and Genomic, but
nevertheless, the precision is significantly higher for all groups except
Locus and Protein (p-value < 0.01 in Welch’s t-test) with respect to empty and
random models (Fig. 3, panel d). This suggests that the model trained with the
gene attributes is more reliable, as it obtains fewer false positives with
similar recall. Thus, for PhenoLinker we incorporate both phenotype and all
gene attributes, as we see that this improves the results.

###  6.5 PhenoLinker predicts valid phenotype-gene associations in an
alternative dataset

To add more evidence of the model’s capabilities beyond using HPO as the
source for validation, we used an additional source of information to further
assess the quality of PhenoLinker predictions: the Gene2Phenotype dataset.

Gene2Phenotype (G2P) [34] is an online system containing evidence-based
datasets linking genes, disease and phenotypes. In these datasets, there are
relationships that have not been included in HPO. Therefore we could use them
to test the efficiency of our model in predicting these new associations and
to see if the model is able to infer them just from the HPO information. G2P
is divided into six disease panels: Cancer, Cardiac, Developmental Disorder
(DD), Eye, Skeletal and Skin. We ruled out the Cancer panel of our analyses as
the size of its associations is rather small. Out of the 5829 associations in
Gene2Phenotype not included in HPO, 3366 are scored by PhenoLinker with
confidence to be included into HPO (i.e., their PhenoLinker score is greater
than 0.5).

To get a proper assessment of whether PhenoLinker behaves better than random
chance on the G2P associations, we determined whether PhenoLinker scores
generated on G2P associations are higher than random ones. For such purpose,
we defined two levels of randomness: (1) full random phenotype-gene pairs when
randomly sampling from the whole HPO set of genes and phenotypes, not
necessarily included into G2P and (2) random phenotype-gene pairs created by
randomly sampling phenotypes and genes included in the corresponding G2P
panel. Supplementary Table S3 illustrates the fold enrichment analysis we
obtain in the comparison. For the full random case, the best fold is obtained
for the skeletal panel. PhenoLinker score is 19 times higher when working on
genuine associations than on fully random ones. The lowest is found in the DD
panel, with a fold of 14. When comparing genuine associations with those
randomly created using phenotypes and genes used by G2P, the highest fold is
of 4.9 for the DD panel, while the lowest is found in the skeletal panel with
a value of 2.5.

After a careful inspection of the PhenoLinker predictions with high confidence
(>= 0.9 score) on the DD panel links, we found really remarkable hits. For
example, our model found that PTEN is associated with delayed speech and
language development (HP:0000750). This relationship is supported by several
publications [41, 42]. In fact, this relationship has been added in a later
version of HPO too. As another interesting example, our model also found that
RPS23 is associated with intellectual disability (HP:0001249). This is also
supported by the study presented in [43].

###  6.6 PhenoLinker predictions improve the efficiency of a genetic
diagnostic application

Thanks to the application of PhenoLinker, we have been able to improve the
results of a genetic diagnostic software that prioritizes variants to assist
geneticists in finding the causal variants within the exome data of clinically
diagnosed patients pending the genetic diagnosis. The software in its basic
form assigns scores to genetic variants based on the probability that they are
causing the phenotypes identified in the patient. By adding our predictions to
the system, specifically those predictions with a score greater than 0.95 and
reanalyzing real cases with pending genetic definitive diagnosis we have
managed to identify a number of additional variants previously validated by
clinicians and genetic experts. In Table 3 we show these 11 gene-phenotype
associations and the literature that supported these relationships.

Gene | New phenotypes | Evidence
---|---|---
TUBGCP2 | Neurodevelopmental delay (HP:0012758), Global developmental delay (HP:0001263) | [44, 45]
POF1B | Hypogonadism (HP:0000135), Hypothyroidism (HP:0000821) | [46]
EXOC7 | Hypotonia (HP:0001252), Delayed speech and language development (HP:0000750), Language impairment (HP:0002463) | [47]
KDM3B | Global developmental delay (HP:0001263) | [48, 49]
TRAPPC10 | Neurodevelopmental delay (HP:0012758), Global developmental delay (HP:0001263), Delayed speech and language development (HP:0000750) | [50, 51]
NRXN1 | Neurodevelopmental delay (HP:0012758), Global developmental delay (HP:0001263) | [52, 53]
MRPS28 | Neurodevelopmental delay (HP:0012758), Global developmental delay (HP:0001263) | [54]
UCP2 | Obesity (HP:0001513) | [55]
OPA1 | Microcephaly (HP:0000252) | [56]
ARHGEF10 | Peripheral neuropathy (HP:0009830) | [57]
ZNF335 | Aplasia/Hypoplasia of the corpus callosum (HP:0007370) | [58]

Table 3: New gene-phenotype associations were found using the PhenoLinker
model. It shows for each gene the new associated phenotypes and their
descriptions, with the literature supporting the relation.

###  6.7 Results and explanations are available online

An interactive online application for PhenoLinker is available to consult the
results of a specific gene set and their corresponding linked phenotypes in
Hugging Face Spaces, in this link. The application is focused on the genes
from the intellectual disability Genomics England amber panel [33], and
phenotypes have been selected such that the gene-phenotype association is not
found in HPO and the score produced by the model is greater than 0.9. We
selected this gene set as an interesting collection of potential genes
associated with intellectual disability. PhenoLinker predicts that out of the
380 genes in the panel, 143 are linked to at least one HPO phenotype term
associated with intellectual disability (Sup. Tables S2 and 4). One of those
143 genes is CACNB4. To illustrate the interpretability provided by our model,
the Web application presents a detailed examination of the gene-phenotype
relationship involving CACNB4 and the term Intellectual Disability
(HP:0001249). PhenoLinker assigned a high score of 0.93 to this association,
and the Web application presents the results using a bar chart that elucidates
the interpretation of the findings based on the attributes of the nodes with
most relevant contribution to the phenotype-gene link (see Fig 4). Notably,
this gene exhibits significant expression levels in the cerebellum, as
supported by our data and GTEx analysis (see Sup. Fig. 1). This expression
level emerged as a crucial feature in scoring this association. Furthermore,
our investigation is substantiated by relevant literature discussing variants
in CACNB4, intellectual disability, and cerebellar atrophy [59], and by the
inclusion of this relationship in a later version of HPO (April 5, 2023).

To complement link interpretation based on contributions from each node
attribute, the Web application includes an interactive 3-D representation of
the embeddings for gene and phenotype, collapsed into three dimensions using
t-SNE. In this way, users can visually explore the phenotype associations of a
selected gene or the gene associations of a selected phenotype via a color map
that highlights the relevant scores. See Fig. 4 for a summary of what is shown
at the Web app for CACNB4.

Figure 4: a) 3D representation by TSNE of the embeddings of the phenotypes,
and by a heat map of the scores of these phenotypes in relation to the CACNB4
gene. It is available in the shiny application, where you can interact with
the plot and select the phenotypes to see their characteristics and the
specific score. b) Example of interpretation obtained for CACNB4 and phenotype
HP:0001249 (Intellectual Disability) in the shiny application. The 10 most
important attributes and their weights are shown in a bar chart.

##  7 Discussion

In this article, we have introduced a heterogeneous structure graph neural
network (HSGNN) designed for predicting gene-phenotype associations using, as
training data, HPO and gene descriptions based on expression, genomics and
protein interaction information, that is called PhenoLinker. The key
characteristics of PhenoLinker can be summarized as follows: (1) Graph Neural
Networks (GNNs): We harness GNNs to model relationships and patterns within
the gene-phenotype graph. GNNs excel in capturing neighborhood information and
propagating it efficiently. (2) Heterogeneous Information Networks (HINs):
This implies that our network consists of nodes belonging to various classes
with distinct attributes. This versatility allows us to simultaneously
consider different types of connections between genes and phenotypes. (3)
Utilization of specific gene attributes: This encompasses genomic and
transcriptomic data sourced from a variety of origins. It enriches our model
with comprehensive and diverse information for each gene. (4) Incorporating
phenotype attributes: We achieve this by employing a pre-trained BioBERT to
transform the detailed phenotype descriptions provided by HPO into a numerical
vector of attributes using the generated embeddings.

Comprehensive experiments, which encompassed repeated random sub-sampling
validation (Fig. 2 and Fig. 3), validation similar to others (Table 1), and
temporal validation (Table 2), lead us to get better predictions than state of
the art methods. However, when we evaluated PhenoLinker´s gene-phenotype
predictions like others did, to allow a fair method comparison, there are
three important differences that should be considered: (1) We utilized genes
instead of proteins, aligning with the HPO’s representation of gene-phenotype
associations. In contrast, others like HPODNets [5] employed proteins and
mapped genes to proteins using the UniProt ID mapping tool, subsequently
assigning identical associations to all proteins derived from a single gene.
This approach is necessary for them due to their use of protein networks, a
distinction from our methodology. (2) The HPO Version: We utilized the
December 15, 2022 edition of HPO, whereas HPODNets employed the March 27,
2020, version of HPO. As of now, the March 27, 2020, version is no longer
accessible. (3) Imbalanced Training: We incorporated significantly fewer false
gene-phenotype associations for training our model compared to others (4 vs 26
negative sampling per positive). This discrepancy could have a detrimental
effect on our evaluation metrics because, when we evaluated our model in a
similar manner to others, we had to assess a considerably larger number of
false associations than expected or that had been previously trained. Despite
this, our results outperformed those of the previous methods.

When we assessed our temporal validation performance, we achieved superior
results compared to previous methods (73.3% for AUCPR in our model vs 11.5%
for M-AUCPR and 8.3% for m-AUCPR in HPODNets using the averages of the results
in both cases). However, it is essential to consider that different HPO
versions were utilized due to the varying dates. We employed 7 versions of HPO
that were released subsequent to the one utilized for training the model (from
2023-01-27 to 2023-10-09), and we assessed the performance of the model with
each of these versions. This process enabled us to gather insights into the
model’s future performance and the connection between the inclusion and
removal of gene-phenotype associations and its performance. Aside from the
initial subsequent HPO version (2023-01-27), we distinctly observe a decrease
in AUCPR over time (0.79 to 0.70). In terms of the impact of including and
removing gene-phenotype associations, it is evident that as the negative ratio
increases, the AUCPR decreases (Table 2). For example, in the HPO version
(2023-01-27), we observe an AUCPR of 0.59 with a negative ratio of 0.45,
whereas in the HPO version (2023-06-17), we see an AUCPR of 0.80 with a
reduced negative ratio of 0.12. Hence, we can deduce that the model performs
better when considering new gene-phenotype associations as opposed to
evaluating previously bad-included, unfavorable gene-phenotype associations in
HPO. This phenomenon may be attributed to the characteristics of the training
process.

The utility of PhenoLinker was further demonstrated by using valid, and
reliable, phenotype-gene associations found at Gene2Phenotype but not
available in HPO yet (Fig. 4). Out of 5829 associations, 3366 should be
considered for inclusion in HPO following PhenoLinker. We also incorporated
our predictions as additional annotations of variants in cases pending of a
genetic disease diagnosis previously validated by clinicians (Table 3). By
incorporating gene-phenotype associations identified by our model into a
production software tool that leverages gene-phenotype information, we were
able to identify an additional 11 causative variants not discovered before and
confirmed by experts. Hence, this suggests that PhenoLinker may be useful for
variant-ranking tools like PhenoVar [60], Phrank [61], Starvar [62], Xrare
[63], DeepPVP [64] and commercial solutions (e.g., Emedgene, or Moon
platforms).

An advantage of our model is its capacity to streamline the analysis of genes
that have not yet been integrated into HPO, thanks to the inclusion of gene
attributes into the model. HPO currently encompasses 4619 genes, whereas there
are approximately 20,000 protein-coding genes in existence [65]. It is highly
probable that some of the unused genes will manifest unknown links with
phenotypes in HPO. At the very least, we believe that it is worth exploring
this possibility in the future and assessing the accuracy of the highest
quality predictions of special interest because of the phenotypes involved,
through wet lab experiments, population analysis and community research.

Another novel aspect introduced in this paper was the utilization of a pre-
trained BioBERT model to incorporate phenotype descriptions as attributes. We
showcased that using these text descriptions coded as embeddings enhances the
model’s performance. But using embeddings poses new challenges. On the one
hand, the BioBERT is not perfectly tailored or fine-tuned for this specific
problem. We will explore fine-tuning to improve the embeddings and hopefully
the PhenoLinker predictions. On the other hand, it would be enormously useful
to be able to assign semantics to the R768 dimensions of the embeddings. This
would considerably enrich the explanatory capabilities of PhenoLinker on its
gene-phenotype associations on top of what we already do with gene attributes
(see Fig. 4, panel b and the Web app). We believe that exploring BioBERT’s
explainability could offer valuable insights and enhance the quality of our
work.

For the training and evaluation of this model, a server with an NVIDIA Tesla
A100 with 40GB of VRAM, an Intel(R) Xeon Silver 4316 CPU @ 2.30GHz and 78 GB
of RAM has been used. On this hardware, the average time of a training epoch
is 30 seconds, model inference is on the order of milliseconds and the
production time of an interpretation is in the order of seconds. The size of
our model is shown in the Supplementary Fig. S2, where we can see the number
of parameters and the layers of the model.

##  8 Conclusions

We have created a heterogeneous structure graph neural network called
PhenoLinker specifically tailored for predicting gene-phenotype associations
using HPO, and it has exhibited superior performance compared to previously
published models. The model’s predictions, scores, explanations, and various
visualization methods are accessible to the research community via a Hugging
Face Space (link). Finally, our model’s gene-phenotype predictions have
contributed to enhanced genetic diagnoses in real-life cases.

## References

  * [1] J. Jumper, R. Evans, A. Pritzel, T. Green, M. Figurnov, O. Ronneberger, K. Tunyasuvunakool, R. Bates, A. Žídek, A. Potapenko, A. Bridgland, C. Meyer, S. A. A. Kohl, A. J. Ballard, A. Cowie, B. Romera-Paredes, S. Nikolov, R. Jain, J. Adler, T. Back, S. Petersen, D. Reiman, E. Clancy, M. Zielinski, M. Steinegger, M. Pacholska, T. Berghammer, S. Bodenstein, D. Silver, O. Vinyals, A. W. Senior, K. Kavukcuoglu, P. Kohli, and D. Hassabis, “Highly accurate protein structure prediction with AlphaFold,” Nature, vol. 596, pp. 583–589, Aug. 2021.
  * [2] OpenAI, “GPT-4 Technical Report,” Mar. 2023.
  * [3] J. Cheng, G. Novati, J. Pan, C. Bycroft, A. Žemgulytė, T. Applebaum, A. Pritzel, L. H. Wong, M. Zielinski, T. Sargeant, R. G. Schneider, A. W. Senior, J. Jumper, D. Hassabis, P. Kohli, and Ž. Avsec, “Accurate proteome-wide missense variant effect prediction with AlphaMissense,” Science, vol. 381, p. eadg7492, Sept. 2023.
  * [4] X.-M. Zhang, L. Liang, L. Liu, and M.-J. Tang, “Graph Neural Networks and Their Current Applications in Bioinformatics,” Frontiers in Genetics, vol. 12, 2021.
  * [5] L. Liu, H. Mamitsuka, and S. Zhu, “HPODNets: Deep graph convolutional networks for predicting human protein–phenotype associations,” Bioinformatics, vol. 38, pp. 799–808, Jan. 2022.
  * [6] Z. Zhang, L. Chen, F. Zhong, D. Wang, J. Jiang, S. Zhang, H. Jiang, M. Zheng, and X. Li, “Graph neural network approaches for drug-target interactions,” Current Opinion in Structural Biology, vol. 73, p. 102327, Apr. 2022.
  * [7] A. R. Deans, S. E. Lewis, E. Huala, S. S. Anzaldo, M. Ashburner, J. P. Balhoff, D. C. Blackburn, J. A. Blake, J. G. Burleigh, B. Chanet, L. D. Cooper, M. Courtot, S. Csösz, H. Cui, W. Dahdul, S. Das, T. A. Dececchi, A. Dettai, R. Diogo, R. E. Druzinsky, M. Dumontier, N. M. Franz, F. Friedrich, G. V. Gkoutos, M. Haendel, L. J. Harmon, T. F. Hayamizu, Y. He, H. M. Hines, N. Ibrahim, L. M. Jackson, P. Jaiswal, C. James-Zorn, S. Köhler, G. Lecointre, H. Lapp, C. J. Lawrence, N. Le Novère, J. G. Lundberg, J. Macklin, A. R. Mast, P. E. Midford, I. Mikó, C. J. Mungall, A. Oellrich, D. Osumi-Sutherland, H. Parkinson, M. J. Ramírez, S. Richter, P. N. Robinson, A. Ruttenberg, K. S. Schulz, E. Segerdell, K. C. Seltmann, M. J. Sharkey, A. D. Smith, B. Smith, C. D. Specht, R. B. Squires, R. W. Thacker, A. Thessen, J. Fernandez-Triana, M. Vihinen, P. D. Vize, L. Vogt, C. E. Wall, R. L. Walls, M. Westerfeld, R. A. Wharton, C. S. Wirkner, J. B. Woolley, M. J. Yoder, A. M. Zorn, and P. Mabee, “Finding Our Way through Phenotypes,” PLoS Biology, vol. 13, p. e1002033, Jan. 2015.
  * [8] D. Botstein and N. Risch, “Discovering genotypes underlying human phenotypes: Past successes for mendelian disease, future approaches for complex disease,” Nature Genetics, vol. 33, pp. 228–237, Mar. 2003.
  * [9] Y. Li and J. C. Patra, “Genome-wide inferring gene–phenotype relationship by walking on the heterogeneous network,” Bioinformatics, vol. 26, pp. 1219–1224, May 2010.
  * [10] A. Brand, H. Brand, and T. Schulte in den Bäumen, “The impact of genetics and genomics on public health,” European Journal of Human Genetics, vol. 16, pp. 5–13, Jan. 2008.
  * [11] B. R. Korf, R. E. Pyeritz, and W. W. Grody, “3 - Nature and Frequency of Genetic Disease,” in Emery and Rimoin’s Principles and Practice of Medical Genetics and Genomics (Seventh Edition) (R. E. Pyeritz, B. R. Korf, and W. W. Grody, eds.), pp. 47–51, Academic Press, Jan. 2019.
  * [12] S. Nguengang Wakap, D. M. Lambert, A. Olry, C. Rodwell, C. Gueydan, V. Lanneau, D. Murphy, Y. Le Cam, and A. Rath, “Estimating cumulative point prevalence of rare diseases: Analysis of the Orphanet database,” European Journal of Human Genetics, vol. 28, pp. 165–173, Feb. 2020.
  * [13] C. A. Deisseroth, J. Birgmeier, E. E. Bodle, J. N. Kohler, D. R. Matalon, Y. Nazarenko, C. A. Genetti, C. A. Brownstein, K. Schmitz-Abe, K. Schoch, H. Cope, R. Signer, J. A. Martinez-Agosto, V. Shashi, A. H. Beggs, M. T. Wheeler, J. A. Bernstein, and G. Bejerano, “ClinPhen extracts and prioritizes patient phenotypes directly from medical records to expedite genetic disease diagnosis,” Genetics in Medicine, vol. 21, pp. 1585–1593, July 2019.
  * [14] S. Marwaha, J. W. Knowles, and E. A. Ashley, “A guide for the diagnosis of rare and undiagnosed disease: Beyond the exome,” Genome Medicine, vol. 14, p. 23, Feb. 2022.
  * [15] T. Zemojtel, S. Köhler, L. Mackenroth, M. Jäger, J. Hecht, P. Krawitz, L. Graul-Neumann, S. Doelken, N. Ehmke, M. Spielmann, N. C. Øien, M. R. Schweiger, U. Krüger, G. Frommer, B. Fischer, U. Kornak, R. Flöttmann, A. Ardeshirdavani, Y. Moreau, S. E. Lewis, M. Haendel, D. Smedley, D. Horn, S. Mundlos, and P. N. Robinson, “Effective diagnosis of genetic disease by computational phenotype analysis of the disease-associated genome,” Science translational medicine, vol. 6, p. 252ra123, Sept. 2014.
  * [16] C. F. Wright, P. Campbell, R. Y. Eberhardt, S. Aitken, D. Perrett, S. Brent, P. Danecek, E. J. Gardner, V. K. Chundru, S. J. Lindsay, K. Andrews, J. Hampstead, J. Kaplanis, K. E. Samocha, A. Middleton, J. Foreman, R. J. Hobson, M. J. Parker, H. C. Martin, D. R. FitzPatrick, M. E. Hurles, and H. V. Firth, “Genomic Diagnosis of Rare Pediatric Disease in the United Kingdom and Ireland,” New England Journal of Medicine, vol. 388, pp. 1559–1571, Apr. 2023.
  * [17] N. D. Rao, J. Kaganovsky, E. A. Malouf, S. Coe, J. Huey, D. Tsinajinne, S. Hassan, K. M. King, S. M. Fullerton, A. T. Chen, and B. H. Shirts, “Diagnostic yield of genetic screening in a diverse, community-ascertained cohort,” Genome Medicine, vol. 15, p. 26, Apr. 2023.
  * [18] S. Köhler, M. Gargano, N. Matentzoglu, L. C. Carmody, D. Lewis-Smith, N. A. Vasilevsky, D. Danis, G. Balagura, G. Baynam, A. M. Brower, T. J. Callahan, C. G. Chute, J. L. Est, P. D. Galer, S. Ganesan, M. Griese, M. Haimel, J. Pazmandi, M. Hanauer, N. L. Harris, M. J. Hartnett, M. Hastreiter, F. Hauck, Y. He, T. Jeske, H. Kearney, G. Kindle, C. Klein, K. Knoflach, R. Krause, D. Lagorce, J. A. McMurry, J. A. Miller, M. C. Munoz-Torres, R. L. Peters, C. K. Rapp, A. M. Rath, S. A. Rind, A. Z. Rosenberg, M. M. Segal, M. G. Seidel, D. Smedley, T. Talmy, Y. Thomas, S. A. Wiafe, J. Xian, Z. Yüksel, I. Helbig, C. J. Mungall, M. A. Haendel, and P. N. Robinson, “The Human Phenotype Ontology in 2021,” Nucleic Acids Research, vol. 49, pp. D1207–D1217, Jan. 2021.
  * [19] 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, vol. 25, pp. 25–29, May 2000.
  * [20] The Gene Ontology Consortium, S. A. Aleksander, J. Balhoff, S. Carbon, J. M. Cherry, H. J. Drabkin, D. Ebert, M. Feuermann, P. Gaudet, N. L. Harris, D. P. Hill, R. Lee, H. Mi, S. Moxon, C. J. Mungall, A. Muruganugan, T. Mushayahama, P. W. Sternberg, P. D. Thomas, K. Van Auken, J. Ramsey, D. A. Siegele, R. L. Chisholm, P. Fey, M. C. Aspromonte, M. V. Nugnes, F. Quaglia, S. Tosatto, M. Giglio, S. Nadendla, G. Antonazzo, H. Attrill, G. dos Santos, S. Marygold, V. Strelets, C. J. Tabone, J. Thurmond, P. Zhou, S. H. Ahmed, P. Asanitthong, D. Luna Buitrago, M. N. Erdol, M. C. Gage, M. Ali Kadhum, K. Y. C. Li, M. Long, A. Michalak, A. Pesala, A. Pritazahra, S. C. C. Saverimuttu, R. Su, K. E. Thurlow, R. C. Lovering, C. Logie, S. Oliferenko, J. Blake, K. Christie, L. Corbani, M. E. Dolan, H. J. Drabkin, D. P. Hill, L. Ni, D. Sitnikov, C. Smith, A. Cuzick, J. Seager, L. Cooper, J. Elser, P. Jaiswal, P. Gupta, P. Jaiswal, S. Naithani, M. Lera-Ramirez, K. Rutherford, V. Wood, J. L. De Pons, M. R. Dwinell, G. T. Hayman, M. L. Kaldunski, A. E. Kwitek, S. J. F. Laulederkind, M. A. Tutaj, M. Vedi, S.-J. Wang, P. D’Eustachio, L. Aimo, K. Axelsen, A. Bridge, N. Hyka-Nouspikel, A. Morgat, S. A. Aleksander, J. M. Cherry, S. R. Engel, K. Karra, S. R. Miyasato, R. S. Nash, M. S. Skrzypek, S. Weng, E. D. Wong, E. Bakker, T. Z. Berardini, L. Reiser, A. Auchincloss, K. Axelsen, G. Argoud-Puy, M.-C. Blatter, E. Boutet, L. Breuza, A. Bridge, C. Casals-Casas, E. Coudert, A. Estreicher, M. Livia Famiglietti, M. Feuermann, A. Gos, N. Gruaz-Gumowski, C. Hulo, N. Hyka-Nouspikel, F. Jungo, P. Le Mercier, D. Lieberherr, P. Masson, A. Morgat, I. Pedruzzi, L. Pourcel, S. Poux, C. Rivoire, S. Sundaram, A. Bateman, E. Bowler-Barnett, H. Bye-A-Jee, P. Denny, A. Ignatchenko, R. Ishtiaq, A. Lock, Y. Lussi, M. Magrane, M. J. Martin, S. Orchard, P. Raposo, E. Speretta, N. Tyagi, K. Warner, R. Zaru, A. D. Diehl, R. Lee, J. Chan, S. Diamantakis, D. Raciti, M. Zarowiecki, M. Fisher, C. James-Zorn, V. Ponferrada, A. Zorn, S. Ramachandran, L. Ruzicka, and M. Westerfield, “The Gene Ontology knowledgebase in 2023,” Genetics, vol. 224, p. iyad031, May 2023.
  * [21] M. Kanehisa and S. Goto, “KEGG: Kyoto Encyclopedia of Genes and Genomes,” Nucleic Acids Research, vol. 28, pp. 27–30, Jan. 2000.
  * [22] L. M. Schriml, J. B. Munro, M. Schor, D. Olley, C. McCracken, V. Felix, J. A. Baron, R. Jackson, S. M. Bello, C. Bearer, R. Lichenstein, K. Bisordi, N. C. Dialo, M. Giglio, and C. Greene, “The Human Disease Ontology 2022 update,” Nucleic Acids Research, vol. 50, pp. D1255–D1261, Jan. 2022.
  * [23] C. Peng, S. Dieck, A. Schmid, A. Ahmad, A. Knaus, M. Wenzel, L. Mehnert, B. Zirn, T. Haack, S. Ossowski, M. Wagner, T. Brunet, N. Ehmke, M. Danyel, S. Rosnev, T. Kamphans, G. Nadav, N. Fleischer, H. Fröhlich, and P. Krawitz, “CADA: Phenotype-driven gene prioritization based on a case-enriched knowledge graph,” NAR Genomics and Bioinformatics, vol. 3, Sept. 2021.
  * [24] M. Notaro, M. Schubach, P. N. Robinson, and G. Valentini, “Prediction of Human Phenotype Ontology terms by means of hierarchical ensemble methods,” BMC Bioinformatics, vol. 18, p. 449, Oct. 2017.
  * [25] R. Patel, Y. Guo, A. Alhudhaif, F. Alenezi, S. A. Althubiti, and K. Polat, “Graph-Based Link Prediction between Human Phenotypes and Genes,” Mathematical Problems in Engineering, vol. 2022, p. e7111647, Mar. 2022.
  * [26] F. Shen, S. Peng, Y. Fan, A. Wen, S. Liu, Y. Wang, L. Wang, and H. Liu, “HPO2Vec+: Leveraging heterogeneous knowledge resources to enrich node embeddings for the Human Phenotype Ontology,” Journal of Biomedical Informatics, vol. 96, p. 103246, Aug. 2019.
  * [27] J. A. Botía, S. Guelfi, D. Zhang, K. D’Sa, R. Reynolds, D. Onah, E. M. McDonagh, A. R. Martin, A. Tucci, A. Rendon, H. Houlden, J. Hardy, and M. Ryten, “G2P: Using machine learning to understand and predict genes causing rare neurological disorders,” Mar. 2018.
  * [28] M. Lek, K. J. Karczewski, E. V. Minikel, K. E. Samocha, E. Banks, T. Fennell, A. H. O’Donnell-Luria, J. S. Ware, A. J. Hill, B. B. Cummings, T. Tukiainen, D. P. Birnbaum, J. A. Kosmicki, L. E. Duncan, K. Estrada, F. Zhao, J. Zou, E. Pierce-Hoffman, J. Berghout, D. N. Cooper, N. Deflaux, M. DePristo, R. Do, J. Flannick, M. Fromer, L. Gauthier, J. Goldstein, N. Gupta, D. Howrigan, A. Kiezun, M. I. Kurki, A. L. Moonshine, P. Natarajan, L. Orozco, G. M. Peloso, R. Poplin, M. A. Rivas, V. Ruano-Rubio, S. A. Rose, D. M. Ruderfer, K. Shakir, P. D. Stenson, C. Stevens, B. P. Thomas, G. Tiao, M. T. Tusie-Luna, B. Weisburd, H.-H. Won, D. Yu, D. M. Altshuler, D. Ardissino, M. Boehnke, J. Danesh, S. Donnelly, R. Elosua, J. C. Florez, S. B. Gabriel, G. Getz, S. J. Glatt, C. M. Hultman, S. Kathiresan, M. Laakso, S. McCarroll, M. I. McCarthy, D. McGovern, R. McPherson, B. M. Neale, A. Palotie, S. M. Purcell, D. Saleheen, J. M. Scharf, P. Sklar, P. F. Sullivan, J. Tuomilehto, M. T. Tsuang, H. C. Watkins, J. G. Wilson, M. J. Daly, D. G. MacArthur, and Exome Aggregation Consortium, “Analysis of protein-coding genetic variation in 60,706 humans,” Nature, vol. 536, pp. 285–291, Aug. 2016.
  * [29] J. Harrow, A. Frankish, J. M. Gonzalez, E. Tapanari, M. Diekhans, F. Kokocinski, B. L. Aken, D. Barrell, A. Zadissa, S. Searle, I. Barnes, A. Bignell, V. Boychenko, T. Hunt, M. Kay, G. Mukherjee, J. Rajan, G. Despacio-Reyes, G. Saunders, C. Steward, R. Harte, M. Lin, C. Howald, A. Tanzer, T. Derrien, J. Chrast, N. Walters, S. Balasubramanian, B. Pei, M. Tress, J. M. Rodriguez, I. Ezkurdia, J. van Baren, M. Brent, D. Haussler, M. Kellis, A. Valencia, A. Reymond, M. Gerstein, R. Guigó, and T. J. Hubbard, “GENCODE: The reference human genome annotation for The ENCODE Project,” Genome Research, vol. 22, pp. 1760–1774, Sept. 2012.
  * [30] A. Busch and K. J. Hertel, “HEXEvent: A database of Human EXon splicing Events,” Nucleic Acids Research, vol. 41, pp. D118–D124, Jan. 2013.
  * [31] D. Szklarczyk, A. Franceschini, S. Wyder, K. Forslund, D. Heller, J. Huerta-Cepas, M. Simonovic, A. Roth, A. Santos, K. P. Tsafou, M. Kuhn, P. Bork, L. J. Jensen, and C. von Mering, “STRING v10: Protein-protein interaction networks, integrated over the tree of life,” Nucleic Acids Research, vol. 43, pp. D447–452, Jan. 2015.
  * [32] P. Deka, A. Jurek-Loughrey, and D. P, “Evidence Extraction to Validate Medical Claims in Fake News Detection,” in Health Information Science (A. Traina, H. Wang, Y. Zhang, S. Siuly, R. Zhou, and L. Chen, eds.), Lecture Notes in Computer Science, (Cham), pp. 3–15, Springer Nature Switzerland, 2022.
  * [33] “Intellectual disability - microarray and sequencing (Version 5.337).” https://panelapp.genomicsengland.co.uk/panels/285/.
  * [34] A. Thormann, M. Halachev, W. McLaren, D. J. Moore, V. Svinti, A. Campbell, S. M. Kerr, M. Tischkowitz, S. E. Hunt, M. G. Dunlop, M. E. Hurles, C. F. Wright, H. V. Firth, F. Cunningham, and D. R. FitzPatrick, “Flexible and scalable diagnostic filtering of genomic variants using G2P with Ensembl VEP,” Nature communications, vol. 10, p. 2373, May 2019.
  * [35] W. L. Hamilton, R. Ying, and J. Leskovec, “Inductive Representation Learning on Large Graphs,” Sept. 2018.
  * [36] A. L. Maas, A. Y. Hannun, and A. Y. Ng, “Rectifier Nonlinearities Improve Neural Network Acoustic Models,” 2013.
  * [37] M. Sundararajan, A. Taly, and Q. Yan, “Axiomatic Attribution for Deep Networks,” in Proceedings of the 34th International Conference on Machine Learning, pp. 3319–3328, PMLR, July 2017.
  * [38] N. Kokhlikyan, V. Miglani, M. Martin, E. Wang, B. Alsallakh, J. Reynolds, A. Melnikov, N. Kliushkina, C. Araya, S. Yan, and O. Reblitz-Richardson, “Captum: A unified and generic model interpretability library for PyTorch,” Sept. 2020.
  * [39] L. van der Maaten and G. Hinton, “Visualizing Data using t-SNE,” Journal of Machine Learning Research, vol. 9, no. 86, pp. 2579–2605, 2008.
  * [40] D. Garcia Gasulla, C. U. Cortés García, E. Ayguadé Parra, and J. J. Labarta Mancho, “Evaluating link prediction on large graphs,” in Artificial Intelligence Research and Development: Proceedings of the 18th International Conference of the Catalan Association for Artificial Intelligence, pp. 90–99, IOS Press, 2015.
  * [41] K. Cummings, A. Watkins, C. Jones, R. Dias, and A. Welham, “Behavioural and psychological features of PTEN mutations: A systematic review of the literature and meta-analysis of the prevalence of autism spectrum disorder characteristics,” Journal of Neurodevelopmental Disorders, vol. 14, p. 1, Jan. 2022.
  * [42] M. Steele, M. Uljarević, G. Rached, T. W. Frazier, J. M. Phillips, R. A. Libove, R. M. Busch, P. Klaas, J. A. Martinez-Agosto, S. Srivastava, C. Eng, M. Sahin, and A. Y. Hardan, “Psychiatric Characteristics Across Individuals With PTEN Mutations,” Frontiers in Psychiatry, vol. 12, p. 672070, Aug. 2021.
  * [43] N. A. Paolini, M. Attwood, S. B. Sondalle, C. M. d. S. Vieira, A. M. van Adrichem, F. M. di Summa, M.-F. O’Donohue, P.-E. Gleizes, S. Rachuri, J. W. Briggs, R. Fischer, P. J. Ratcliffe, M. W. Wlodarski, R. H. Houtkooper, M. von Lindern, T. W. Kuijpers, J. D. Dinman, S. J. Baserga, M. E. Cockman, and A. W. MacInnes, “A Ribosomopathy Reveals Decoding Defective Ribosomes Driving Human Dysmorphism,” The American Journal of Human Genetics, vol. 100, pp. 506–522, Mar. 2017.
  * [44] S. Gungor, Y. Oktay, S. Hiz, Á. Aranguren-Ibáñez, I. Kalafatcilar, A. Yaramis, E. Karaca, U. Yis, E. Sonmezler, B. Ekinci, M. Aslan, E. Yilmaz, B. Özgör, S. Balaraju, N. Szabo, S. Laurie, S. Beltran, D. G. MacArthur, D. Hathazi, A. Töpf, A. Roos, H. Lochmuller, I. Vernos, and R. Horvath, “Autosomal recessive variants in TUBGCP2 alter the $\gamma$-tubulin ring complex leading to neurodevelopmental disease,” iScience, vol. 24, p. 101948, Jan. 2021.
  * [45] T. Mitani, J. Punetha, I. Akalin, D. Pehlivan, M. Dawidziuk, Z. Coban Akdemir, S. Yilmaz, E. Aslan, J. V. Hunter, H. Hijazi, C. M. Grochowski, S. N. Jhangiani, E. Karaca, J. M. Fatih, P. Iwanowski, T. Gambin, P. Wlasienko, A. Goszczanska-Ciuchta, M. Bekiesinska-Figatowska, M. Hosseini, S. Arzhangi, H. Najmabadi, J. A. Rosenfeld, H. Du, D. Marafi, S. Blaser, R. Teitelbaum, R. Silver, J. E. Posey, H.-H. Ropers, R. A. Gibbs, W. Wiszniewski, J. R. Lupski, D. Chitayat, K. Kahrizi, and P. Gawlinski, “Bi-allelic Pathogenic Variants in TUBGCP2 Cause Microcephaly and Lissencephaly Spectrum Disorders,” American Journal of Human Genetics, vol. 105, pp. 1005–1015, Nov. 2019.
  * [46] A. Lacombe, H. Lee, L. Zahed, M. Choucair, J.-M. Muller, S. F. Nelson, W. Salameh, and E. Vilain, “Disruption of POF1B Binding to Nonmuscle Actin Filaments Is Associated with Premature Ovarian Failure,” The American Journal of Human Genetics, vol. 79, pp. 113–119, July 2006.
  * [47] M. E. Coulter, D. Musaev, E. M. DeGennaro, X. Zhang, K. Henke, K. N. James, R. S. Smith, R. S. Hill, J. N. Partlow, Muna Al-Saffar, A. S. Kamumbu, N. Hatem, A. J. Barkovich, J. Aziza, N. Chassaing, M. S. Zaki, T. Sultan, L. Burglen, A. Rajab, L. Al-Gazali, G. H. Mochida, M. P. Harris, J. G. Gleeson, and C. A. Walsh, “Regulation of human cerebral cortical development by EXOC7 and EXOC8, components of the exocyst complex, and roles in neural progenitor cell proliferation and survival,” Genetics in Medicine, vol. 22, pp. 1040–1050, June 2020.
  * [48] I. J. Diets, R. van der Donk, K. Baltrunaite, E. Waanders, M. R. F. Reijnders, A. J. M. Dingemans, R. Pfundt, A. T. Vulto-van Silfhout, L. Wiel, C. Gilissen, J. Thevenon, L. Perrin, A. Afenjar, C. Nava, B. Keren, S. Bartz, B. Peri, G. Beunders, N. Verbeek, K. van Gassen, I. Thiffault, M. Cadieux-Dion, L. Huerta-Saenz, M. Wagner, V. Konstantopoulou, J. Vodopiutz, M. Griese, A. Boel, B. Callewaert, H. G. Brunner, T. Kleefstra, N. Hoogerbrugge, B. B. A. de Vries, V. Hwa, A. Dauber, J. Y. Hehir-Kwa, R. P. Kuiper, and M. C. J. Jongmans, “De Novo and Inherited Pathogenic Variants in KDM3B Cause Intellectual Disability, Short Stature, and Facial Dysmorphism,” The American Journal of Human Genetics, vol. 104, pp. 758–766, Apr. 2019.
  * [49] M. Tabaku, S. Tomori, P. Cullufi, E. Dervishi, O. Paknia, and P. Bauer, “A novel de novo pathogenic variant in KDM3B gene at the first Albanian case of Diets-Jongmans syndrome: DIJOS,” Molecular Genetics and Metabolism Reports, vol. 33, p. 100927, Oct. 2022.
  * [50] L. E. Rawlins, H. Almousa, S. Khan, S. C. Collins, M. P. Milev, J. Leslie, D. Saint-Dic, V. Khan, A. M. Hincapie, J. O. Day, L. McGavin, C. Rowley, G. V. Harlalka, V. E. Vancollie, W. Ahmad, C. J. Lelliott, A. Gul, B. Yalcin, A. H. Crosby, M. Sacher, and E. L. Baple, “Biallelic variants in TRAPPC10 cause a microcephalic TRAPPopathy disorder in humans and mice,” PLoS Genetics, vol. 18, p. e1010114, Mar. 2022.
  * [51] M. Sacher, N. Shahrzad, H. Kamel, and M. P. Milev, “TRAPPopathies: An emerging set of disorders linked to variations in the genes encoding transport protein particle (TRAPP)-associated proteins,” Traffic, vol. 20, no. 1, pp. 5–26, 2019.
  * [52] M. S. L. Ching, Y. Shen, W.-H. Tan, S. S. Jeste, E. M. Morrow, X. Chen, N. M. Mukaddes, S.-Y. Yoo, E. Hanson, R. Hundley, C. Austin, R. E. Becker, G. T. Berry, K. Driscoll, E. C. Engle, S. Friedman, J. F. Gusella, F. M. Hisama, M. B. Irons, T. Lafiosca, E. LeClair, D. T. Miller, M. Neessen, J. D. Picker, L. Rappaport, C. M. Rooney, D. P. Sarco, J. M. Stoler, C. A. Walsh, R. R. Wolff, T. Zhang, R. H. Nasir, B.-L. Wu, and Children’s Hospital Boston Genotype Phenotype Study Group, “Deletions of NRXN1 (neurexin-1) predispose to a wide spectrum of developmental disorders,” American Journal of Medical Genetics. Part B, Neuropsychiatric Genetics: The Official Publication of the International Society of Psychiatric Genetics, vol. 153B, pp. 937–947, June 2010.
  * [53] M. Sciacca, L. Marino, G. Vitaliti, R. Falsaperla, and S. Marino, “NRXN1 Deletion in Two Twins’ Genotype and Phenotype: A Clinical Case and Literature Review,” Children, vol. 9, p. 698, May 2022.
  * [54] J. Pulman, B. Ruzzenente, L. Bianchi, M. Rio, N. Boddaert, A. Munnich, A. Rötig, and M. D. Metodiev, “Mutations in the MRPS28 gene encoding the small mitoribosomal subunit protein bS1m in a patient with intrauterine growth retardation, craniofacial dysmorphism and multisystemic involvement,” Human Molecular Genetics, vol. 28, pp. 1445–1462, May 2019.
  * [55] J. Li, R. Jiang, X. Cong, and Y. Zhao, “UCP2 gene polymorphisms in obesity and diabetes, and the role of UCP2 in cancer,” FEBS letters, vol. 593, pp. 2525–2534, Sept. 2019.
  * [56] A. Nasca, T. Rizza, M. Doimo, A. Legati, A. Ciolfi, D. Diodato, C. Calderan, G. Carrara, E. Lamantea, C. Aiello, M. Di Nottia, M. Niceta, C. Lamperti, A. Ardissone, S. Bianchi-Marzoli, G. Iarossi, E. Bertini, I. Moroni, M. Tartaglia, L. Salviati, R. Carrozzo, and D. Ghezzi, “Not only dominant, not only optic atrophy: Expanding the clinical spectrum associated with OPA1 mutations,” Orphanet Journal of Rare Diseases, vol. 12, p. 89, May 2017.
  * [57] G. K. Boora, A. A. Kulkarni, R. Kanwar, P. Beyerlein, R. Qin, M. S. Banck, K. J. Ruddy, J. Pleticha, C. A. Lynch, R. J. Behrens, S. Züchner, C. L. Loprinzi, and A. S. Beutler, “Association of the Charcot-Marie-Tooth disease gene ARHGEF10 with paclitaxel induced peripheral neuropathy in NCCTG N08CA (Alliance),” Journal of the Neurological Sciences, vol. 357, pp. 35–40, Oct. 2015.
  * [58] K. Stouffs, A. B. Stergachis, T. Vanderhasselt, A. Dica, S. Janssens, L. Vandervore, A. Gheldof, O. Bodamer, K. Keymolen, S. Seneca, I. Liebaers, D. Jayaraman, H. E. Hill, J. N. Partlow, C. A. Walsh, and A. C. Jansen, “Expanding the clinical spectrum of biallelic ZNF335 variants,” Clinical genetics, vol. 94, pp. 246–251, Aug. 2018.
  * [59] P. Coste de Bagneaux, L. von Elsner, T. Bierhals, M. Campiglio, J. Johannsen, G. J. Obermair, M. Hempel, B. E. Flucher, and K. Kutsche, “A homozygous missense variant in CACNB4 encoding the auxiliary calcium channel beta4 subunit causes a severe neurodevelopmental disorder and impairs channel and non-channel functions,” PLoS Genetics, vol. 16, p. e1008625, Mar. 2020.
  * [60] Y. J. Trakadis, C. Buote, J.-F. Therriault, P.-É. Jacques, H. Larochelle, and S. Lévesque, “PhenoVar: A phenotype-driven approach in clinical genomics for the diagnosis of polymalformative syndromes,” BMC Medical Genomics, vol. 7, p. 22, May 2014.
  * [61] K. A. Jagadeesh, J. Birgmeier, H. Guturu, C. A. Deisseroth, A. M. Wenger, J. A. Bernstein, and G. Bejerano, “Phrank measures phenotype sets similarity to greatly improve Mendelian diagnostic disease prioritization,” Genetics in Medicine, vol. 21, pp. 464–470, Feb. 2019.
  * [62] S. Kafkas, M. Abdelhakim, M. Uludag, A. Althagafi, M. Alghamdi, and R. Hoehndorf, “Starvar: Symptom-based tool for automatic ranking of variants using evidence from literature and genomes,” BMC Bioinformatics, vol. 24, p. 294, July 2023.
  * [63] Q. Li, K. Zhao, C. D. Bustamante, X. Ma, and W. H. Wong, “Xrare: A machine learning method jointly modeling phenotypes and genetic evidence for rare disease diagnosis,” Genetics in Medicine: Official Journal of the American College of Medical Genetics, vol. 21, pp. 2126–2134, Sept. 2019.
  * [64] I. Boudellioua, M. Kulmanov, P. N. Schofield, G. V. Gkoutos, and R. Hoehndorf, “DeepPVP: Phenotype-based prioritization of causative variants using deep learning,” BMC Bioinformatics, vol. 20, p. 65, Feb. 2019.
  * [65] S. L. Salzberg, “Open questions: How many genes do we have?,” BMC Biology, vol. 16, p. 94, Aug. 2018.

Supplementary Materials

Group | Source | Attributes | Description
---|---|---|---
|  | ExACpLI | The probability of being loss-of-function, lof, intolerant of both heterozygous and homozygous lof variants
|  | ExACpRec | The probability of being intolerant of homozygous, but not heterozygous lof variants
|  | ExACpNull | The probability of being tolerant of both heterozygous and homozygous lof variants
ExAC | The Exome Aggregation Consortium database | ExACpMiss | Corrected missense Z score, taking into account that higher Z scores indicate that the transcript is more intolerant of variation
|  | CountsOverlap | The number of genes overlapping a gene of interest
Locus | R package GRanges | CountsProtCodOverlap | The number of overlapping protein coding genes
Protein | STRING database | StringCombined | The number of catalogued interactions for the corresponding proteins
|  | constitutiveexons | The constituve exons of each gene
|  | ESTcount | Number of detected ESTs
|  | alt3EST | 3’ Alternative EST
|  | alt5EST | 5’ Alternative EST
|  | alt3.5EST | Alternative 3’ and 5’ EST
|  | GeneLength | Gene length
|  | TranscriptCount | Number of transcripts for the corresponding gene
|  | GCcontent | Percentage of GC content
|  | NumJunctions | The total number of unique exon-exon junctions generated by each gene
Genomic |  "Ensembl version 72" | IntronicLength | The intronic length of each gene
Expression | GTEx V6 gene expression dataset | 47 tissue expression attributes | Boolean that indicates if the gene is expressed in the tissue or not

Table S1: Attribute groups for genes, where we can see the sources for each group, the attributes and a description for each one. Phenotype description | HPO ID
---|---
Cognitive impairment | HP:0100543
Mental deterioration | HP:0001268
Intellectual disability | HP:0001249
Global developmental delay | HP:0001263
Intellectual disability, severe | HP:0010864
Intellectual disability, progressive | HP:0006887
Intellectual disability, profound | HP:0002187
Intellectual disability, moderate | HP:0002342
Intellectual disability, mild | HP:0001256
Intellectual disability, borderline | HP:0006889
Table S2: List of HPO phenotype terms associated with intellectual disability.

Disease Panel | Relations | Best Threshold Disease Related | Best Fold Enrichment Disease Related | Best Threshold Full Randomized | Best Fold Enrichment Full Randomized
---|---|---|---|---|---
Developmental Disorder (DD) | 3030 | 0,9 | 4,9 | 0,85 | 14
Skeletal | 436 | 0,95 | 2,5 | 0,85 | 19
Eye | 1470 | 0,93 | 3,7 | 0,93 | 15,1
Skin | 767 | 0,92 | 3 | 0,9 | 16,3

Table S3: Fold enrichment analysis. For each panel, the maximum fold reached and the threshold where it is reached are shown. This analysis is performed both by comparing with disease-related random distributions and with completely random distributions. Figure S1: Bulk tissue gene expression for gene CACNB4, where it is shown that the gene is mainly expressed in the cerebellar hemisphere. Figure S2: Summary of the PhenoLinker model. The layers, input/output shapes and number of parameters are shown for an inference on the test set with a batch of 2048 links. Table S4: Gene-phenotype associations, where genes are in the intellectual disability Genomics England amber panel and phenotypes are selected from Supplementary Table S1. Only associations with a score >0.9 are shown. Gene | Phenotype | Name | Score
---|---|---|---
ACADSB | HP:0001249 | Intellectual disability | 0,9821
ACADVL | HP:0001249 | Intellectual disability | 0,9540
ACADVL | HP:0001263 | Global developmental delay | 0,9672
ACAT1 | HP:0100543 | Cognitive impairment | 0,9037
ACTA2 | HP:0001263 | Global developmental delay | 0,9839
ACVR1 | HP:0001263 | Global developmental delay | 0,9285
ADAMTS10 | HP:0001263 | Global developmental delay | 0,9875
AIMP2 | HP:0010864 | Intellectual disability, severe | 0,9079
AIMP2 | HP:0100543 | Cognitive impairment | 0,9314
ALDOA | HP:0001263 | Global developmental delay | 0,9156
ALPL | HP:0001249 | Intellectual disability | 0,9369
ALPL | HP:0001263 | Global developmental delay | 0,9626
ALX3 | HP:0001263 | Global developmental delay | 0,9944
ANTXR1 | HP:0001263 | Global developmental delay | 0,9864
ATP11A | HP:0001249 | Intellectual disability | 0,9785
B4GALT1 | HP:0001249 | Intellectual disability | 0,9902
BSND | HP:0001263 | Global developmental delay | 0,9942
CACNA2D2 | HP:0001249 | Intellectual disability | 0,9930
CACNA2D2 | HP:0001268 | Mental deterioration | 0,9043
CACNA2D2 | HP:0100543 | Cognitive impairment | 0,9125
CACNB4 | HP:0001249 | Intellectual disability | 0,9306
CACNB4 | HP:0001263 | Global developmental delay | 0,9497
CACNB4 | HP:0100543 | Cognitive impairment | 0,9673
CASR | HP:0001249 | Intellectual disability | 0,9455
CASR | HP:0001263 | Global developmental delay | 0,9408
CCDC88A | HP:0010864 | Intellectual disability, severe | 0,9540
CCDC88A | HP:0100543 | Cognitive impairment | 0,9559
CDKN1C | HP:0001249 | Intellectual disability | 0,9861
CLCN2 | HP:0001249 | Intellectual disability | 0,9090
CLCN2 | HP:0001263 | Global developmental delay | 0,9191
CNPY3 | HP:0100543 | Cognitive impairment | 0,9339
COPB1 | HP:0001263 | Global developmental delay | 0,9943
COPB1 | HP:0100543 | Cognitive impairment | 0,9015
COQ9 | HP:0001249 | Intellectual disability | 0,9870
CPSF3 | HP:0010864 | Intellectual disability, severe | 0,9237
CPSF3 | HP:0100543 | Cognitive impairment | 0,9328
CTC1 | HP:0001249 | Intellectual disability | 0,9734
CYP27A1 | HP:0001263 | Global developmental delay | 0,9893
CYP2U1 | HP:0001268 | Mental deterioration | 0,9022
DALRD3 | HP:0001256 | Intellectual disability, mild | 0,9088
DALRD3 | HP:0010864 | Intellectual disability, severe | 0,9366
DDOST | HP:0001249 | Intellectual disability | 0,9617
DENND5A | HP:0001249 | Intellectual disability | 0,9921
DENND5A | HP:0010864 | Intellectual disability, severe | 0,9354
DENND5A | HP:0100543 | Cognitive impairment | 0,9431
DLG1 | HP:0001263 | Global developmental delay | 0,9879
DOHH | HP:0001268 | Mental deterioration | 0,9284
DPH2 | HP:0001249 | Intellectual disability | 0,9820
DPM3 | HP:0001263 | Global developmental delay | 0,9366
DYNC1I2 | HP:0100543 | Cognitive impairment | 0,9342
EMG1 | HP:0001249 | Intellectual disability | 0,9681
ENTPD1 | HP:0001263 | Global developmental delay | 0,9935
EXOC2 | HP:0001249 | Intellectual disability | 0,9951
EXOC2 | HP:0010864 | Intellectual disability, severe | 0,9307
EXOC7 | HP:0001249 | Intellectual disability | 0,9834
EXOC7 | HP:0010864 | Intellectual disability, severe | 0,9051
EXOSC8 | HP:0001249 | Intellectual disability | 0,9948
EXOSC8 | HP:0001268 | Mental deterioration | 0,9011
EXOSC8 | HP:0100543 | Cognitive impairment | 0,9544
FARSB | HP:0001249 | Intellectual disability | 0,9761
FDFT1 | HP:0001249 | Intellectual disability | 0,9879
FGF13 | HP:0010864 | Intellectual disability, severe | 0,9321
FGF14 | HP:0001263 | Global developmental delay | 0,9622
FGF14 | HP:0001268 | Mental deterioration | 0,9503
FRAS1 | HP:0001263 | Global developmental delay | 0,9954
FREM2 | HP:0001263 | Global developmental delay | 0,9929
FRMD5 | HP:0001268 | Mental deterioration | 0,9123
FRMD5 | HP:0100543 | Cognitive impairment | 0,9756
FTO | HP:0001249 | Intellectual disability | 0,9910
GBA2 | HP:0001263 | Global developmental delay | 0,9940
GCSH | HP:0001263 | Global developmental delay | 0,9815
GCSH | HP:0100543 | Cognitive impairment | 0,9410
GEMIN4 | HP:0001249 | Intellectual disability | 0,9805
GJB1 | HP:0001249 | Intellectual disability | 0,9662
GJB1 | HP:0001263 | Global developmental delay | 0,9620
GJB3 | HP:0001263 | Global developmental delay | 0,9622
GLS | HP:0001249 | Intellectual disability | 0,9912
GLS | HP:0001268 | Mental deterioration | 0,9296
GLS | HP:0100543 | Cognitive impairment | 0,9644
GOT2 | HP:0001263 | Global developmental delay | 0,9952
GPSM2 | HP:0001263 | Global developmental delay | 0,9795
GRIA1 | HP:0100543 | Cognitive impairment | 0,9540
GRM7 | HP:0001256 | Intellectual disability, mild | 0,9193
GRM7 | HP:0002342 | Intellectual disability, moderate | 0,9058
GRM7 | HP:0010864 | Intellectual disability, severe | 0,9439
GRM7 | HP:0100543 | Cognitive impairment | 0,9739
GSS | HP:0001263 | Global developmental delay | 0,9312
GSX2 | HP:0001249 | Intellectual disability | 0,9912
GSX2 | HP:0001268 | Mental deterioration | 0,9305
GSX2 | HP:0100543 | Cognitive impairment | 0,9125
HADHB | HP:0001249 | Intellectual disability | 0,9670
HEATR3 | HP:0001263 | Global developmental delay | 0,9729
IFT27 | HP:0001263 | Global developmental delay | 0,9359
IFT43 | HP:0001249 | Intellectual disability | 0,9597
IFT43 | HP:0001263 | Global developmental delay | 0,9659
IREB2 | HP:0001249 | Intellectual disability | 0,9896
IREB2 | HP:0100543 | Cognitive impairment | 0,9469
ISCA2 | HP:0001249 | Intellectual disability | 0,9904
ISCA2 | HP:0001263 | Global developmental delay | 0,9962
ISCA2 | HP:0100543 | Cognitive impairment | 0,0946
ITGA7 | HP:0001263 | Global developmental delay | 0,9977
KCNC3 | HP:0001268 | Mental deterioration | 0,9514
KIF4A | HP:0001263 | Global developmental delay | 0,9729
LARS2 | HP:0001249 | Intellectual disability | 0,9686
LARS2 | HP:0001263 | Global developmental delay | 0,9817
LIPT2 | HP:0001249 | Intellectual disability | 0,9920
LIPT2 | HP:0001268 | Mental deterioration | 0,9391
LIPT2 | HP:0100543 | Cognitive impairment | 0,9033
LNPK | HP:0001249 | Intellectual disability | 0,0985
MAN2C1 | HP:0001263 | Global developmental delay | 0,9922
MAPKAPK5 | HP:0001249 | Intellectual disability | 0,9872
MED12L | HP:0001263 | Global developmental delay | 0,9983
MPV17 | HP:0001249 | Intellectual disability | 0,9597
NAGS | HP:0001249 | Intellectual disability | 0,9397
NBAS | HP:0001249 | Intellectual disability | 0,9720
NBAS | HP:0001263 | Global developmental delay | 0,9835
NBN | HP:0001263 | Global developmental delay | 0,9857
NCAPD2 | HP:0001263 | Global developmental delay | 0,0984
NCAPG2 | HP:0001249 | Intellectual disability | 0,9921
NDUFAF1 | HP:0001249 | Intellectual disability | 0,9929
NDUFAF1 | HP:0100543 | Cognitive impairment | 0,9037
NDUFAF2 | HP:0001268 | Mental deterioration | 0,9626
NDUFAF2 | HP:0100543 | Cognitive impairment | 0,0982
NDUFAF5 | HP:0001268 | Mental deterioration | 0,0962
NDUFAF5 | HP:0100543 | Cognitive impairment | 0,9754
NECAP1 | HP:0001256 | Intellectual disability, mild | 0,9380
NECAP1 | HP:0010864 | Intellectual disability, severe | 0,0931
NFIB | HP:0001249 | Intellectual disability | 0,9537
NFIB | HP:0001263 | Global developmental delay | 0,9796
NPHP3 | HP:0001249 | Intellectual disability | 0,9645
NUP188 | HP:0001249 | Intellectual disability | 0,9862
NUP188 | HP:0001263 | Global developmental delay | 0,9964
NUP214 | HP:0001249 | Intellectual disability | 0,0959
NUP62 | HP:0001263 | Global developmental delay | 0,9959
NUP62 | HP:0001268 | Mental deterioration | 0,9730
NUP62 | HP:0100543 | Cognitive impairment | 0,9743
PAM16 | HP:0001249 | Intellectual disability | 0,9423
PDE10A | HP:0001263 | Global developmental delay | 0,9713
PIK3C2A | HP:0001249 | Intellectual disability | 0,9808
PLXNA1 | HP:0001249 | Intellectual disability | 0,9885
PNPO | HP:0001249 | Intellectual disability | 0,9914
POMK | HP:0010864 | Intellectual disability, severe | 0,9448
PPFIBP1 | HP:0010864 | Intellectual disability, severe | 0,9259
PRKACB | HP:0001263 | Global developmental delay | 0,9870
PRKAR1B | HP:0001249 | Intellectual disability | 0,9935
PRKAR1B | HP:0001256 | Intellectual disability, mild | 0,9154
PRKD1 | HP:0001249 | Intellectual disability | 0,9791
PRRT2 | HP:0001263 | Global developmental delay | 0,9863
PRRT2 | HP:0001268 | Mental deterioration | 0,9842
PRRT2 | HP:0100543 | Cognitive impairment | 0,9937
PSMB8 | HP:0001263 | Global developmental delay | 0,9835
PTH1R | HP:0001249 | Intellectual disability | 0,9678
PTH1R | HP:0001263 | Global developmental delay | 0,9826
RBPJ | HP:0001263 | Global developmental delay | 0,9931
RNF220 | HP:0001263 | Global developmental delay | 0,9923
RNF220 | HP:0001268 | Mental deterioration | 0,9214
RNF220 | HP:0100543 | Cognitive impairment | 0,9016
RPS23 | HP:0001249 | Intellectual disability | 0,9773
RPS23 | HP:0001263 | Global developmental delay | 0,9942
RSPRY1 | HP:0001263 | Global developmental delay | 0,9886
RUSC2 | HP:0001256 | Intellectual disability, mild | 0,9198
SACS | HP:0001263 | Global developmental delay | 0,9944
SACS | HP:0001268 | Mental deterioration | 0,9445
SACS | HP:0100543 | Cognitive impairment | 0,9435
SCN1B | HP:0001268 | Mental deterioration | 0,0092
SCN1B | HP:0010864 | Intellectual disability, severe | 0,9118
SEC31A | HP:0001249 | Intellectual disability | 0,9899
SEC31A | HP:0001263 | Global developmental delay | 0,9978
SEMA6B | HP:0001268 | Mental deterioration | 0,9262
SEMA6B | HP:0100543 | Cognitive impairment | 0,9407
SHROOM4 | HP:0001263 | Global developmental delay | 0,9927
SLC26A2 | HP:0001249 | Intellectual disability | 0,9675
SLC26A2 | HP:0001263 | Global developmental delay | 0,9840
SLC35D1 | HP:0001263 | Global developmental delay | 0,9137
SLC39A13 | HP:0001249 | Intellectual disability | 0,9155
SLC39A13 | HP:0001263 | Global developmental delay | 0,9344
SLC5A5 | HP:0001263 | Global developmental delay | 0,9794
SLC5A7 | HP:0001263 | Global developmental delay | 0,9995
SLC9A7 | HP:0001249 | Intellectual disability | 0,9569
SMAD3 | HP:0001249 | Intellectual disability | 0,9838
SMAD3 | HP:0001263 | Global developmental delay | 0,9969
SMG9 | HP:0010864 | Intellectual disability, severe | 0,9150
SUCLA2 | HP:0100543 | Cognitive impairment | 0,9392
TAB2 | HP:0001263 | Global developmental delay | 0,9975
TBX1 | HP:0100543 | Cognitive impairment | 0,9128
TGFB1 | HP:0001249 | Intellectual disability | 0,9711
THRB | HP:0001249 | Intellectual disability | 0,9328
THRB | HP:0001263 | Global developmental delay | 0,9265
TKFC | HP:0001249 | Intellectual disability | 0,9879
TMEM147 | HP:0001263 | Global developmental delay | 0,9985
TNIK | HP:0001263 | Global developmental delay | 0,9879
TNIK | HP:0100543 | Cognitive impairment | 0,9497
TNR | HP:0001249 | Intellectual disability | 0,9923
TNR | HP:0100543 | Cognitive impairment | 0,0915
TRAPPC10 | HP:0002342 | Intellectual disability, moderate | 0,9173
TRAPPC11 | HP:0100543 | Cognitive impairment | 0,9324
TRAPPC2L | HP:0001249 | Intellectual disability | 0,9785
TSEN15 | HP:0001268 | Mental deterioration | 0,9669
TSEN15 | HP:0002187 | Intellectual disability, profound | 0,9638
TSEN15 | HP:0006887 | Intellectual disability, progressive | 0,9357
TSEN15 | HP:0010864 | Intellectual disability, severe | 0,9618
TSEN15 | HP:0100543 | Cognitive impairment | 0,9809
TSPOAP1 | HP:0001263 | Global developmental delay | 0,9527
TSPOAP1 | HP:0001268 | Mental deterioration | 0,9634
TSPOAP1 | HP:0100543 | Cognitive impairment | 0,9796
TUBGCP2 | HP:0001249 | Intellectual disability | 0,9929
TUBGCP2 | HP:0001263 | Global developmental delay | 0,9986
TUBGCP2 | HP:0010864 | Intellectual disability, severe | 0,9300
TUBGCP2 | HP:0100543 | Cognitive impairment | 0,9392
UFC1 | HP:0001263 | Global developmental delay | 0,9983
VIPAS39 | HP:0001249 | Intellectual disability | 0,9661
VPS33B | HP:0001249 | Intellectual disability | 0,9715
VPS50 | HP:0001249 | Intellectual disability | 0,9929
VPS51 | HP:0001249 | Intellectual disability | 0,9937
VPS51 | HP:0010864 | Intellectual disability, severe | 0,9133
WASHC5 | HP:0100543 | Cognitive impairment | 0,9016
XPA | HP:0001263 | Global developmental delay | 0,9820