Integration of Human Protein Sequence and Protein-Protein Interaction Data by Graph Autoencoder to Identify Novel Protein-Abnormal Phenotype Associations

Abstract

Understanding gene functions and their associated abnormal phenotypes is crucial in the prevention, diagnosis and treatment against diseases. The Human Phenotype Ontology (HPO) is a standardized vocabulary for describing the phenotype abnormalities associated with human diseases. However, the current HPO annotations are far from completion, and only a small fraction of human protein-coding genes has HPO annotations. Thus, it is necessary to predict protein-phenotype associations using computational methods. Protein sequences can indicate the structure and function of the proteins, and interacting proteins are more likely to have same function. It is promising to integrate these features for predicting HPO annotations of human protein. We developed GraphPheno, a semi-supervised method based on graph autoencoders, which does not require feature engineering to capture deep features from protein sequences, while also taking into account the topological properties in the protein-protein interaction network to predict the relationships between human genes/proteins and abnormal phenotypes. Cross validation and independent dataset tests show that GraphPheno has satisfactory prediction performance. The algorithm is further confirmed on automatic HPO annotation for no-knowledge proteins under the benchmark of the second Critical Assessment of Functional Annotation, 2013-2014 (CAFA2), where GraphPheno surpasses most existing methods. Further bioinformatics analysis shows that predicted certain phenotype-associated genes using GraphPheno share similar biological properties with known ones. In a case study on the phenotype of abnormality of mitochondrial respiratory chain, top prioritized genes are validated by recent papers. We believe that GraphPheno will help to reveal more associations between genes and phenotypes, and contribute to the discovery of drug targets.

Keywords: deep learning; graph autoencoder; protein-phenotype associations prediction.

Figure 1

Overview of GraphPheno. The model consists of three modules: (A) the Protein embedding module: this module takes proteins interaction and sequence information as input. PPI network was converted to the format of adjacency matrix A. Proteins amino acid sequences were embedded using the conjoint triad (CT) method and served as proteins initial features X; (B) the VGAE module for feature embedding: This module consists of a two-layer GCN encoder and a dot product decoder, and generates latent representations Z based on both topological information from PPI network and protein sequence features. The adjacency matrix $$\tilde{A}$$ is reconstructed using the latent variable Z through the dot product decoder. (C) Neural network module for prediction: Gold standard dataset is used to train this module. This module takes VGAE embedding as an input to produce prediction scores for each gene-phenotype association.

Figure 2

Performance evaluation of GraphPheno. (A,B) Performance comparison of various prediction models in term-centric evaluation (A) and protein-centric evaluation (B). Confidence intervals (95%) were determined using bootstrapping with 100 iterations. (C) The distributions of AUCs for the prediction of 3741 HPO terms by GraphPheno model against five-fold cross-validation prediction. (D) The distributions of AUCs for the prediction of 2993 HPO terms by GraphPheno model against independent test set validation. The mean value of AUCs are plotted in dotted lines.

Figure 3

Performance comparison under the benchmark of the CAFA2 challenge using F-max. GraphPheno (green) was compared with the top performing CAFA2 participating methods (light gray), baseline methods in CAFA2 challenge (red for Naive, and blue for BLAST), and several HPO predicting methods proposed after CAFA2 challenge (dark gray). F-max is the maximum value of F-measure over all thresholds. Confidence intervals (95%) were determined using bootstrapping with 100 iterations.

Figure 4

Predicted and known annotated genes of 4369 phenotypes share similar biological properties. Box plots of the sequence consistency (A), number of protein–protein interactions (B), gene expression correlation coefficient (C), number of proteins in the smallest shared biological process (D) between predicted genes and random genes with known annotated genes for each phenotype. The random genes were randomly selected from unannotated genes of each phenotype with an equal number of predicted genes. (In the box plots, the middle bar represents the median, and the box represents the interquartile range; bars extend to 1.5× the interquartile range. p-values are calculated by the Student’s t-test and shown on the top of the boxes. *** p-value < 0.001, **** p-value < 0.0001).

Figure 5

Functional analysis of predicted and known annotated genes with the phenotype of Decreased activity of mitochondrial respiratory chain. (A) Overlap between predicted and known annotated “Decreased activity of mitochondrial respiratory chain”-associated genes with respect to enriched GO terms and KEGG pathways. A relatively large shared GO enrichment terms and KEGG enrichment pathways was found. (*** p-value < 0.001, hypergeometric test). Predicted “Decreased activity of mitochondrial respiratory chain”-associated genes are enriched in mitochondrial related GO terms (B) and Mitochondrial related diseases such as Alzheimer’s disease, Parkinson’s disease, etc. (C).

Figure 6

Both predicted and known annotated “Decreased activity of mitochondrial respiratory chain”-associated genes tend to be significantly down regulated in multiple neurodegenerative diseases. Enrichment ratio was calculated as the GeneRatio divided by Background Ratio. GeneRatio refers to the number of predicted or known annotated “Decreased activity of mitochondrial respiratory chain”-associated genes which are significantly down regulated in the GEO dataset divided by the total number of predicted or known annotated “Decreased activity of mitochondrial respiratory chain”-associated genes. Background Ratio refers to the number of significantly down regulated genes in the GEO dataset divided by the total number of genes identified in the GEO dataset. (ns: p-value > 0.05, *** p-value < 0.001, hypergeometric test).

Figure 7

Biological insight into the predicted “Decreased activity of mitochondrial respiratory chain”-associated genes. The oxidative phosphorylation (OXPHOS) system is embedded in the lipid bilayer of the inner mitochondrial membrane (IMM) and is composed of five protein enzyme complexes and two mobile electron carriers namely ubiquinone (CoQ) and cytochrome c (Cyt C). Translocator of the outer and inner mitochondrial membrane (TOM and TIM, respectively) were also shown. Predicted “Decreased activity of mitochondrial respiratory chain”-associated genes were presented in dotted boxes, in which circle and rectangle denote genes which function as the subunits and assemble factors of mitochondrial respiratory chain complex I–V, TOM40 complex and TIM23 complex, respectively. Predicted “Decreased activity of mitochondrial respiratory chain”-related genes validated by recent papers or the newest version of HPO database were marked with asterisks.