Unveiling new disease, pathway, and gene associations via multi-scale neural network

Diseases involve complex modifications to the cellular machinery. The gene expression profile of the affected cells contains characteristic patterns linked to a disease. Hence, new biological knowledge about a disease can be extracted from these profiles, improving our ability to diagnose and assess disease risks. This knowledge can be used for drug re-purposing, or by physicians to evaluate a patient’s condition and co-morbidity risk. Here, we consider differential gene expressions obtained by microarray technology for patients diagnosed with various diseases. Based on these data and cellular multi-scale organization, we aim at uncovering disease–disease, disease–gene and disease–pathway associations. We propose a neural network with structure based on the multi-scale organization of proteins in a cell into biological pathways. We show that this model is able to correctly predict the diagnosis for the majority of patients. Through the analysis of the trained model, we predict disease–disease, disease–pathway, and disease–gene associations and validate the predictions by comparisons to known interactions and literature search, proposing putative explanations for the predictions.


Introduction
A disease is often described by symptoms and affected tissues.However, to give a definite diagnosis, physicians often need to analyse patient samples (e.g., blood samples, or biopsies) for characteristic disease indicators, commonly referred to as disease biomarkers.These may include disregulated genes, or pathways [16,15].Taking into consideration the history of a patient's past and present conditions identifying genetic predispositions, as well as considering associations between diseases, aid in achieving accurate diagnostics and treatments [27].By also taking advantage of the increasing availability of large scale molecular data, precision medicine aims at improving the understanding of the molecular base of diseases on individual basis, as well as the relationships between different conditions [19,33].The benefits from such work are multiple and include drug re-purposing and identification of new disease biomarkers to improve treatments and diagnoses.
Many studies have investigated disease-gene and disease-pathway associations with the objective of improving diagnoses [2,24,54,26].For instance, Zhao et al. [54] propose a ranking of disease genes based on gene expression and protein interactions using Katz-centrality.Hong et al. [26] design a tool that identifies significantly disrupted pathways by comparing patient gene expression against controls collected from other experiments.Cogswell et al. [9] identify putative gene and pathway biomarkers through change in miRNA in Alzheimer's disease.In specific cancers, Abeel et al. [2] use support vector machines and ensemble feature selection methods to select putative gene biomarkers.Ciucci et al. [8] developed a general purpose algorithm based on the analysis of PCA loadings that can be used to identify genes that discriminate between conditions from expression data.
A key issue is that most of these studies consider diseases in isolation, i.e. comparing patients having a disease of interest to healthy individuals while the predicted biomarkers could be shared between various diseases.This limits the discriminative potential of such studies for accurate diagnoses.Indeed, network medicine has shown that diseases can share significant molecular background, as evidenced by numerous studies based on patient historical records [25,27,28], biological knowledge of the diseases [19,33,22], or patient gene expression profiles [44].For instance, Goh et al. [19] build a disease network, which connects diseases that share at least one gene which when mutated is linked to both conditions.Lee et al. [33] construct a disease network of metabolic diseases, connecting pairs of diseases for which associated mutated enzymes catalyze adjacent metabolic reactions.Hidalgo et al. [25] take a different approach by building a disease network based on disease co-morbidities, i.e. two diseases are connected if they tend to co-occur significantly in the patient populations.They used Medicare records of elderly patients to build the network.He et al. [22] propose PCID (Predicting Comorbidity by Integrating Data), an approach to predict disease co-morbidities by aggregating disease similarity scores derived from different data including protein-protein interactions (PPIs), pathways, and functional annotations.
Sánchez-Valle et al. [44] define a disease network, named the Disease Molecular Similarity Network (DMSN) based on patient's differential expression profiles.In their study, the DMSN is generated using positive and negative relative molecular similarities (RR) to measure disease similarity and dissimilarity, respectively, that is then interpreted as an estimate of risk.First, a patient-patient similarity network is generated based on the similarities of patient's differential expression profiles.Next, using the relative similarity score, diseases are related to each other.The resulting network is directed and each edge is associated to a positive or negative label indicating either an increased or decreased risk of developing the target disease if the patient has the source disease.The underlying assumption is that having a given disease can increase the risk of developing a disease characterized by a similar gene expression profile.
In these various approaches, a key issue is that either a single data source is used, such as disease-gene mutational data [19], or no new biological knowledge about a specific disease could be derived from the results (e.g., PCID [22]).
In this work, we propose an integrative framework based on artificial neural networks (NN) to predict diseasedisease links, as well as disease-pathway and disease-gene associations.We train the model to predict patients' diagnoses based on differential gene expression.The NN's structure is designed to mimic the cellular multi-scale functional organization by integrating gene-pathway annotations.This approach follows on from a body of work aiming to build neural networks based on prior information defining the structure of the network [39].For instance, Ma et al. [34] build a neural network using the Gene Ontology [4] directed acyclic graph as a template for the connections.The neural network, named DCell, is then trained to predict phenotype related to cellular fitness from genotype data.The trained DCell predicts cellular growth almost as accurately as laboratory observations.
We show that our framework achieves good predicting performances on our dataset.By analysing the trained NN, i.e. the underlying weight matrices, we show that we can extract biological knowledge relevant for each disease.Specifically, we use the trained NN to predict novel disease-pathway and disease-gene links and from those predictions we extract disease similarity score used to identify putative co-morbidities.We show that our predictions are biologically relevant against established ground-truths and verify the top predictions through manual literature curation ensuring that the sources do not use the same data, to mitigate the risk of argument circularity.

Datasets
We use the original dataset of patient gene expressions used by [44] in which each patient is associated to a single disease (see Supplement 3).For each patient, we define a differential gene expression vector of size corresponding to the number of genes and in which the i th entry is equal to 1, −1, or 0 depending on whether the i th gene is significantly (with 5% cutoff) over-, under-, or normally expressed, respectively, for that patient (see Supplement 3 for details).
Pathway annotations were collected from Reactome database [17] (accessed December 2018).Only the lowest pathways in the hierarchy are considered to avoid dealing with pathway interactions (i.e.pathways containing other pathways).Of those pathways, only the ones that have a Traceable Author Statement (TAS) are kept.In total, we consider 1, 708 pathway annotations.
The final dataset contains 4, 788 samples (patients) diagnosed by one of 83 diseases (see Table 7).In total, 20, 525 genes have their expressions measured, but only a subset is used as input to our method described in the following section, as we restrict ourselves to genes associated to at least one pathway, which leaves 9, 247 genes.

Neural network based data-integration framework
We propose a neural network (NN) predicting a patient diagnosis based on differential gene expression.The structure of the neural network is based on molecular organization, more specifically gene-pathway annotations downloaded from Reactome (see Figure 1).We integrate molecular organization data into our model to reflect the idea that complex diseases, such as cancer, can be the results of the perturbations of groups of genes, as opposed to a single gene.Using Reactome data allows us to incorporate prior knowledge into our model in the form of biologically meaningful groupings of genes.

Phenotypes (Diseases)
Figure 1: Example of neural network architecture.For the first layer, the connections are defined by biological information, i.e. a unit representing a gene is connected to all the biological pathways that the gene is involved in.We do not add any prior knowledge on the last layer, thus it is fully connected.
A neural network can be expressed as a series of matrix multiplications interleaved with non-linear function (See Supplement 3 for more details).Here, we use the softmax function [20] as the last non-linear function of the NN (common choice for multiclass classification problems) and the hyperbolic tangent non-linear function for hidden layers to allow a hidden unit to have values lying in [−1, 1] to represent up-and down-regulations.We use the classical cross-entropy function to define the loss function.Our neural network architecture has only one hidden layer capturing gene-pathway links.Hereafter, we refer to this model as GPD (for Gene-Pathway-Disease).
As reference model, we use a multiclass logistic regression (MLR; no hidden layer).MLR and GPD have a different number of parameters (or free weights): (573, 314) for MLR and (137, 838) for GPD.Note that, due to this imbalance, we do not expect GPD to outperform MLR in the diagnosis prediction task.
To train both GPD and MLR, we first perform a 10-fold cross-validation to fix the number of training epochs.To fix the number of training epochs, we compute the average number of epochs at which the test loss is the smallest across the runs (see Figure 3).As the dataset is imbalanced, we use stratification to split the data, ensuring that at least one patient per disease is in the test set.Using this number of epochs, we perform another 10-fold crossvalidation to evaluate the performance of our models.We use the Adam optimizer [29] with learning rate 0.01 and the layer weights are initialized to small values using the initialization process proposed by He et al. [23].We investigated the addition of classical regularisation techniques -L1, L2, and dropout regularisation -as a mean to reduce the capacity of the model to overfit.We found that the performances were better without any regularisations (see Table 8).
The neural networks are implemented with Tensorflow [1].

Predicting disease-disease, disease-pathway, and disease-gene relationships
To identify associations between diseases and genes or pathways, we perform sensitivity analysis [55] (see Supplement 3).The association score between disease d and unit u (representing a gene, or a pathway) is measured by the intensity of the local variation of the output unit associated with d with respect to perturbation of u.Intuitively, this score measures how the prediction score of disease d is affected by disregulation of a gene/pathway: we quantify the change in one disease score induced by a disregulation in the gene expression, or pathway activation.We test this scoring approach for the prediction of disease-gene and disease-pathway associations.In particular, we rank disease-gene and disease-pathway pairs based on this score and test if the score correlates with known associations through a Precision-Recall and Receiver Operating Characteristic (ROC) analysis.Score thresholds could be inferred from the precision-recall curves.These thresholds would be uninformative in general since they are tied to specific runs and datasets.Hence, we focus on manual validation of the top scoring associations.
Based on this score, we represent each disease by a set formed by the k genes highest scoring genes and a set containing the k pathways highest scoring pathways.We then score disease-disease associations using the Jaccard Index of their sets.The Jaccard Index of two sets S 1 and S 2 is defined as |S1∩S2| |S1∪S2| , where | • | represents the cardinality of a set.Following on from similar approach used in DisGeNET [42], we interpret those associations as co-morbidities.The number of highest scoring pathways and genes considered is set to 150 and 300, respectively, as those numbers gave the best results.

Classification performances
To validate the relevance of our model, we verify that the classification performances are at least on par with competing methods: MLR, Random Forest (RF), Bernoulli Naive Bayes (nB), and Support Vector Machine (SVM) algorithms (we use the implementation available through the scikit-learning python package [41]).We perform 10fold cross-validation for the algorithms to fix the hyperparameters (numbers of trees 100, smoothing parameter 0.001 and penalty parameter 100, respectively) and retain the best performing models in terms of cross-entropy loss (the objective function of the neural networks).
We evaluate performances by computing 3 different scores: cross-entropy loss (CEL), micro-and macro-average precision (Pre µ and Pre M ).Pre µ give a measure of the overall precision of each classifier, while Pre M gives an average of the precisions across the different classes (diseases).The details of each score can be found in Supplement 3. We observe that the neural networks; MLR and GPD, give better, or at least on-par, performances when compared to RF, nB, and SVM classifiers as measured by our three metrics (see Table 1).This observation justifies the relevance of our GPD model.The best model appears to be the multinomial logistic regression (MLR), which corresponds to the most complex neural network model in terms of the number of parameters (or degree of freedom), since MLR has ∼ 4 times more parameters than GPD.This analysis shows that using biological knowledge to guide the structure of neural networks, in the limit of the models proposed, does not improve classification performance compared to the multiclass logistic regression (MLR) and only offers slight improvement when compared to a RF classifier (see Table 1).However, we show, in the following Sections, that the trained GPD models can be more successfully exploited than MLR to extract biological information.Note as well that the gene-pathway information on which GPD relies is both noisy and incomplete, as biological data often is, and that performances should improve as knowledge improves.

Algorithm
Hereafter, we consider for each model (GPD and MLR) the trained NN that gave the lowest cross-entropy loss during the 10-fold cross validation.

Our GPD model uncovers molecular mechanisms of diseases
To uncover molecular mechanisms of disease, i.e., genes and pathways that are associated to specific diseases, we extract predictions from MLR and GPD using the approach described in Section 1.3.We test the performance of our disease-pathway and disease-gene associations predictions by comparing against established databases.We investigate the top predictions of the GPD model through manual search of the literature.

Predicting disease-gene associations
For each model, we compute disease-gene association scores as described in Section1.3and we test the validity of our predictions against DisGeNET database [42].We compare the entire set of predictions against two baselines: the Frequency of Differential Expression (FDE) and the approach introduced by Zhao et al. [54] for de novo disease-gene association prediction (Katz).Those methods are detailed in the Supplement 3.
We use precision-recall and ROC curves to evaluate the performance of our approach and compute the areas under the curves (see Table 2 and Figure 4).Interestingly, we observe that the FDE score is a poor predictor of disease-gene associations.We further observe that GPD is the best performing models for this task with Katz coming second.The relatively poor overall performances can be partially attributed to the incompleteness of the reported disease-gene associations in DisGeNET.To corroborate this hypothesis, we search the literature for support for the top 10 predicted disease-gene associations by the best performing model, GPD (see  We are able to find literature support for 70% of the top 10 predicted disease-gene associations (see Table 3).Furthermore, we find indications that some of our top-scoring, non-validated predictions could be relevant, such as the associations of asthma with UBB and amyotrophic lateral sclerosis (ALS) with PSMD13.Ubuquitin B (UBB) belongs to the ubiquitin-proteasome (UPS) and it is known that aberration in the UPS is responsible for inflammatory and autoimmune diseases such as asthma [53].Moreover, ALS onsets occur typically after age 50 and manifest partially through muscle weakness.PSMD13 is linked to aging [6] and high expression of the gene has been found in skeletal muscle of athletes [31], suggesting that under-expression could be a sign of muscle weakness.
These results validate the relevance of our framework for de novo disease-gene association prediction and confirm the incompleteness of DisGeNET.

Predicting disease-pathway associations
For our GPD model, we compute disease-pathway association scores as described in Section 1.3 and we test the validity of our predictions by comparison with CTD database [13].As a baseline, we consider disease-pathway scores corresponding to the average FDE (AFDE) of genes within the pathway for patients having the disease.
We evaluate the results as done previously for disease-gene associations (see Table 4 and Figure 5).We observe that GPD convincingly outperforms AFDE.The seemingly poor performances of both approaches can partially be attributed to the incompleteness of CTD database.To test this hypothesis, we search the literature for support for the top 10 disease-pathway associations predicted with our GPD (see  We find literature support for 7 out of the top 10 predicted disease-pathway associations (see Table 5).Furthermore, we find indications that some of our top-scoring, non-validated predictions could be relevant, such as the association of autistic disorder with the lactose synthesis pathway (R-HSA-5653890) and the association of Schizophrenia with pathway R-HSA-5683371 linked to microphtalmia.The lactose synthesis pathway (R-HSA-5653890) contains three genes: LALBA, SLC2A1, and B4GALT1.All of those genes might be associated with autistic disorders.One patented method to detect autistic disorder (US20140349977A1) includes LALBA as one of the genes of interest.SLC2A1 mutation has been reported in patients diagnosed with autism [21].And finally, B4GALT1 has been linked with developmental disorders [52].The pathway R-HSA-5683371 is linked to the eye disease microphthalmia.It is known that schizophrenia is linked to eye abnormalities [51].Among the 28 genes involved in that pathway, 12 have been linked to the disease in the literature (GOT2, PDHA1, DLD, GCSH, DLAT, PDHB, DAO, OGDH, DHTKD1, GNMT, DDO, PRODH2).
These results show the relevance of our framework for de novo disease-pathway associations prediction despite relatively low retrieval scores against the ground-truth.

Our GPD model predicts disease-disease relationships
We rank disease-disease pairs based on the score described in Section 1.3 and test our results against a high confidence co-morbidity disease network obtained from a large cohort study by Hidalgo et al. [25].We compare our method against DMSN network [44], restricted to our set of diseases, and three alternatives baselines.For the first alternative, we compute disease-disease association score using our approach defined in Section 1.3 on the trained MLR network, representing each disease by the top 300 highest scoring genes (which gave the best results based on grid search).The final two baselines associate to each disease-disease pair a Jaccard Index score based on 1) the set of genes associated to each disease in DisGeNET [42] and 2) the set of pathways associated to each disease in CTD database [13].The results of the comparison are presented using a precision-recall curve (see Figure 2).
We observe that our approach outperform convincingly the other approaches in the task of retrieving existing comorbidity links between diseases with over 10% increase compared to DMSN and 30% improvement over DisGeNET in terms of area under the precision-recall curve (auprc).These results strongly support our methodology.The scoring based on disease-gene is performing better than disease-pathway, hence we investigate the top 10 scoring disease-disease associations derived from it (see Table 6).
We present and discuss below literature support for the predicted associations between the diseases.Atrial fibrillation has been linked in the literature to thyroid disease [47] which is known to be co-morbid with vitiligo [11].Atrial fibrillation and peripheral vascular disease are well known co-morbid conditions [43].Alcoholism has been linked to the onset of some cancers notably implicating the transcription factor Nanog which itself has been linked to osteosarcoma [35].Additionally, a drug used to treat alcoholism, Disulfiram, has recently been proposed as a potential treatment for osteosarcoma [10].Those information put together suggest shared molecular background   6: Top 10 disease-disease links predicted using our approach based on the trained GPD.
for the two conditions.Rhabdoid cancer is a rare form of aggressive cancer affecting young children and with very poor prognostic, which makes it difficult to evaluate co-morbid conditions.However, rhabdoid cancer is frequently mistaken for medulloblastoma indicating some similarity [7].We found no evidence in the literature supporting a connection between the rare Cornelia de Lange syndrome and Vitiligo.Some studies have observed significant comorbidity between vitiligo and psoriasis and the combination of the two has been linked to cardiovascular diseases, which include peripheral vascular disease [3,11].Atrial fibrillation and cardiac complications have been observed as the result of osteosarcoma [49,14].Leishmaniasis and alcoholic hepatitis are an unlikely co-morbid connection since it would require a patient both to have been infected by parasites of the Leishmania type and have had excessive alcohol intake.However, both disease affects the liver and leishmaniasis has sometimes been misdiagnosed for cirrhosis [18] which suggests that the two diseases might share some similar molecular processes that we would be capturing here.
A case of co-occurence of Sotos syndrome and vitiligo has been reported in the medical literature [12].It has been postulated that non-Hodgkin's lymphoma (which include follicular lymphoma) and osteosarcoma share underlying mechanisms [36].Additionally, miR-202 has been identified as a potential tumor suppressor for both conditions [50].Through this analysis, we have shown that most predicted pairs have either been observed co-occurring or can be connected through underlying mechanisms, thus validating our approach.

Conclusions
In this study, we propose a multi-scale neural-network based framework that integrates gene expression data associated to diseases with gene-pathway information.Our integrative framework allows for simultaneously uncovering novel disease-disease associations and disease molecular mechanisms from patient gene expression profiles through the analysis of trained neural networks.We show that GPD achieve good diagnosis prediction on our dataset showing the validity of our integrative process.Furthermore, we show that the associations predicted from the trained models are biologically meaningful and supported by the literature, thus validating our approach.
While our disease molecular mechanisms are supported by the current knowledge about these diseases, a next step would be to identify among the predicted genes and pathways suitable biomarkers and drug targets that could be used to improve diagnosis, prognosis, and treatment.We leave this for future work.Also, while our multi-scale NN framework integrates the hierarchical functional organization of a cell (from genes to biological pathways), our methodology can be extended to include any dataset pertaining to diseases of interest, e.g., uncovering molecular mechanisms of cancer from patient somatic mutation profiles or linking diseases to non-coding RNA.
Finally, while we focus on patient data with application to diseases, our methodology can be extended to integrate additional omics data with the objective of getting more biologically accurate models for the analyses of patients, tissues, and cells.Some further applications include studies of diseases linked to a specific tissue, studies of cell's specialization, and any study that can benefit from the integration of the hierarchical functional organization of cells.
where Ŷ is the objective, or ground truth, and g(•) is a predefined function.Multinomial logistic regression (MLR) and the proposed GDP architecture can be written as where s is the softmax function, typically used for multiclass classification problems and tanh denotes the hyperbolic tangent.Matrices X and Ŷ represent our data.Each column of X corresponds to the differential gene expression of a patient, and each column of Ŷ corresponds to a patient's diagnosis (the prediction of which is the objective of the framework).W 1 ∈ R n d ×ng and W 2 2 ∈ R n d ×np correspond to fully-connected layers.The layer corresponding to W 2 1 ∈ R nc×ng represents biological pathway membership of the genes, i.e. the trainable weights of the matrix correspond only to entries (i, j) where gene j is part of the i th protein complex.SI Methods Formally, the local variations δf of a single-argument function f due to a change δx = x − x 0 in input can be approximated with the first order Taylor expansion as Thus, the magnitude of the local variations of f with respect to perturbation δx from x 0 is given by | df dx (x 0 )|.Based on this approximation, we extract from each neural network a score between an entity represented by a unit of the neural network (e.g, a pathway, or a gene) and each disease (output unit).Specifically, for a neural network NN, we denote nn i : [0, 1] ni → R no the function corresponding to the operation of a neural network NN from the n i outputs of layer i to the final n o logits of the neural network, i.e. scores before softmax.E.g., for GPD, we have nn 1  2 (x) = W 2 2 tanh W 2 1 x .Then, the association score s i,j,k between the j th output unit of layer i, denoted u i j , of neural network NN, and disease k is given by where the reference point is chosen as the null vector, x 0 = 0, which corresponds implicitly to a healthy state in our formulation.

SI Metrics
The cross-entropy loss (CLE) of a classifier is defined as where m represents the number of samples (patients), n the number of classes (diseases), y ij indicates if patient i is diagnosed with disease j (1 if true 0 otherwise), and ŷij is the j th output value of the classifier for patient i.
A relatively small CEL means that the output probability distribution of a classifier is closer to the deterministic one-hot encoding of the true labeling, i.e. the classifier gives a higher probability to the true class and a very small probability to the other classes.The micro-averaged precision (Pre µ ) of a classifier gives a measure of the overall precision of the classifier and is defined as Pre where tp corresponds to the number of accurately classified patients and m represents the number of patients.The macro-averaged precision (Pre M ) of a classifer gives an average of the precision across the different classes (diseases) and is defined as where tp i corresponds to the number of accurately classified patients for disease i and m i represents the number of patients diagnosed with the same disease.Pre M can be more informative than Pre µ when considering the problem with class imbalance.

SI Baselines
The Frequency of Differential Expression (FDE) score of a disease-gene association corresponds to how frequently that gene is consistently differentially expressed in patients having the disease, i.e., for disease d and gene g, the association score, s dg , is given by where we amalgamate entities (disease, gene, and patient) with their indices, P d denotes the set of patients having disease d, and X corresponds to the data matrix introduced in Methods.The Katz method uses disease specific Protein-Protein Interaction (PPI) network, where each node of a standard PPI network is associated to a score (here the FDE of each gene for the disease considered).The authors then use Katz-centrality on each disease PPI network to extract a final score for each disease-gene association (here we use the absolute value).The higher the score, the higher the association is expected to be true.We download the PPI data from BioGRID [40] and IID [32] and create our PPI network from the union of both databases restricted to our set of genes.Finally, we perform a grid-search to identify the best parameters for the model by trying to maximise the area under the precision-recall curve metric.

Figure 2 :
Figure 2: Precision-recall (top) and ROC (bottom) curves of the test against the disease co-morbidity network built by Hidalgo et al. [25].

Figure 3 :
Figure 3: Train and test loss curve for the MLR model (left) and GDP model (right) with respect to the number of epochs during the cross-validation.The vertical black line indicates the number of epochs which give the lowest loss.

Figure 4 :
Figure 4: Precision-recall curve and ROC curve of our predictions for disease-genes associations.

Figure 5 :
Figure 5: Precision-recall curve and ROC curve of our disease-pathways associations predictions.

Table 1 :
Performances of different classifiers in terms of cross-entropy loss (CEL), micro-and macro-average precisions (Pre µ and Pre M , respectively).Each score is computed across the 10-fold cross-validation and we provide the standard deviation.Bold scores highlight the best scores for each metric.

Table 5 )
. Note that none of those associations are reported in DisGeNET.

Table 2 :
Performance in terms of area under the ROC (AUROC) and area under the precision-recall (AUPRE) for the prediction of disease-gene associations for each methods)

Table 5 )
. Note that none of these associations predicted are reported in CTD database.

Table 4 :
Performance in terms of area under the ROC (AUROC) and area under the precision-recall (AUPRE) for the prediction of disease-pathway associations for each methods)

Table 8 :
Results of cross-validation to fix regularisation hyperparameters (L1-, L2-, or dropout regularisations).The scores correspond to cross-entropy loss.The best results are obtained with no regularisation (score in bold).