HPOAnnotator: improving large-scale prediction of HPO annotations by low-rank approximation with HPO semantic similarities and multiple PPI networks

Background As a standardized vocabulary of phenotypic abnormalities associated with human diseases, the Human Phenotype Ontology (HPO) has been widely used by researchers to annotate phenotypes of genes/proteins. For saving the cost and time spent on experiments, many computational approaches have been proposed. They are able to alleviate the problem to some extent, but their performances are still far from satisfactory. Method For inferring large-scale protein-phenotype associations, we propose HPOAnnotator that incorporates multiple Protein-Protein Interaction (PPI) information and the hierarchical structure of HPO. Specifically, we use a dual graph to regularize Non-negative Matrix Factorization (NMF) in a way that the information from different sources can be seamlessly integrated. In essence, HPOAnnotator solves the sparsity problem of a protein-phenotype association matrix by using a low-rank approximation. Results By combining the hierarchical structure of HPO and co-annotations of proteins, our model can well capture the HPO semantic similarities. Moreover, graph Laplacian regularizations are imposed in the latent space so as to utilize multiple PPI networks. The performance of HPOAnnotator has been validated under cross-validation and independent test. Experimental results have shown that HPOAnnotator outperforms the competing methods significantly. Conclusions Through extensive comparisons with the state-of-the-art methods, we conclude that the proposed HPOAnnotator is able to achieve the superior performance as a result of using a low-rank approximation with a graph regularization. It is promising in that our approach can be considered as a starting point to study more efficient matrix factorization-based algorithms.


Background
Phenotypes refer to observable physical or biological traits of an organism. Revealing the relationships between genes/proteins and their related phenotypes is one of the main objectives of genetics in the post-genome era [1][2][3]. The Human Phenotype Ontology (HPO) [4] is a standardized vocabulary for describing the phenotypic abnormalities associated with human diseases [5]. Being initially populated by using databases of human genes and genetic disorders such as OMIM [6], Orphanet [7] and DECIPHER [8], HPO was later expanded by using literature curation [9]. At present, only small quantities of human protein-coding genes (∼ 3500) have HPO annotations. It is, however, believed that a large number of currently unannotated genes/proteins are related to disease phenotypes. Therefore, it is critical to predict genes/protein-HPO associations by using accurate computational methods.
Currently, HPO contains four sub-ontologies: Organ abnormality, Mode of inheritance, Clinical modifier, and Mortality/Aging. As the main sub-ontology, Organ abnormality describes clinical abnormalities whose first-level children are formed by terms like abnormality of a skeletal system. The Mode of inheritance describes inheritance patterns of phenotypes and contains terms such as Autosomal dominant. The Clinical modifier contains classes that describe typical modifiers of clinical symptoms such as those triggered by carbohydrate ingestion. For Mortality/Aging, it describes the age of death by terms like Death in childhood and Sudden death. The Organ abnormality, Mode of inheritance, Clinical modifier, and Mortality/Aging have ∼12000, 28, 100, and 8 terms, respectively.
The annotations between genes/proteins and HPO terms are very sparse. Specifically, 284621 annotations are for 3459 proteins and 6407 HPO terms with the sparsity of 1.2%. Meanwhile, the annotation growth by time, for example, is about 5%, with adding only 14820 annotations as new ones between June 2017 to December 2017. Since genes/proteins are annotated with multiple HPO terms, the prediction can be regarded as a problem of multi-label predictions. Differing from this, HPO terms, however, form a hierarchical structure. This implies that once a gene/protein is labeled with one HPO term, it should also be labeled with all of its ancestors of this particular HPO term. In other words, when a gene/protein is not labeled with an HPO term, it should not be labeled with all of its descendants, either. That is, general terms are located at the top of the HPO structure, with the term specificity increasing from the root to the leaves. Figure 1 shows a real example of an HPO hierarchical structure (i.e., Directed Acyclic Graph, DAG) and the scale of sub-ontologies.
The existing computational approaches for HPO annotation prediction can be divided into two categories, namely feature-based and network-based methods. The feature-based approaches use gene/protein information as the features to predict its annotations for a query gene/protein. For sparse and noisy data, the incorporation of auxiliary information into original input data generally helps to improve predictive performance. One of these methods, learning to rank, has been demonstrated the superior performance in GO annotation prediction [10], for example. Compared with GO annotations, HPO annotations are, however, more reliable and stable. In addition, the sparseness of HPO annotations is much less than that of GO annotations, with focusing on human proteins and terms under Organ abnormality only. Nevertheless, few existing feature-based models take into consideration HPO information, e.g., the hierarchical structure and cooccurrence of HPO terms. The network-based approaches are more prevalent at present. Usually, multiple networks are integrated into a new large-scale network in order to improve the prediction in these approaches such as random-walk [11] and weighted score computation [12]. However, network-based approaches cannot perform well for sparse data. This is because of disconnected nodes that are commonly encountered in real-world graphs, particulary for sparse data, even though they can be related to each other.
Prediction of the annotations between genes/proteins and HPO terms can be grouped into two categorises: 1) pair prediction, which predicts the missing HPO annotations of existing proteins, and 2) prediction of new proteins, which annotates HPO terms to the totally unannotated proteins. Most existing work belong to the latter category, but few are for the former. To narrow this gap, we focus on the first category in this paper, which is also a famous task in the CAFA challenge. Existing methods for the first category have four major limitations. First, the hierarchy of HPO is completely ignored. The hierarchical structure poses a formidable challenge to a prediction: a model needs to evaluate the associations between a protein and all of its related phenotypes from the deeper levels to the root in the hierarchy. Second, the existing methods do not make full use of the potentials of Protein-Protein Interaction (PPI) networks. For example, a PPI network is modeled in the original annotation space in their models, which may not extract the information effectively. Moreover, multiple PPI networks may be derived from different sources, resulting in the data fusion. Third, only a few known associations are available for training. So they are extremely unbalanced. Specifically, more than half of the terms in HPO are used to annotate zero or only one protein. As a result, such a drastic sparsity makes prediction more challenging. Finally, existing methods usually study the sub-ontologies independently without considering the co-annotations of HPO terms. However, co-annotations are quite common in annotations. It is likely that they help improve prediction results.
To address the above four problems, we apply matrix factorization to approximate a protein-HPO annotation matrix by two factorized low-rank matrices. As such, the latent factors that underlie the HPO annotations can be well captured. Since the HPO annotation matrix is binary, we choose to use Non-negative Matrix Factorization (NMF). NMF has proved to be effective for sparse problems in the field of bioinformatics [13][14][15][16]. Based on our above observations, we propose an NMF-based framework called HPOAnnotator by which to predict missing protein-HPO annotations. In essence, the key idea of our model is to factorize the HPO annotation matrix into two non-negative low-rank latent matrices, which correspond to the respective latent feature spaces of proteins and HPO terms. In addition, the graph Laplacian on PPI networks is performed to exploit their intrinsic geometric structure. Co-annotations and the hierarchical structure of HPO are also incorporated to measure HPO semantic relationships.
We have experimentally validated the performance of HPOAnnotator by comparing it with the three networkbased approaches, which will be reviewed in the related work. The proposed model was tested on the latest largescale HPO data with around 300000 annotations. Experimental results clearly demonstrated that HPOAnnotator outperformed the competing methods under two scenarios: cross-validation and independent test. It indicates that a low-rank approximation and network information are effective for pair prediction. Furthermore, our case studies further provide evidence for the practical use of HPOAnnotator. Note that, the work presented in this paper is the extension of our previous work AiProAnnotator [17] (AiPA for short). The main difference between the two methods is that HPOAnnotator can seamlessly combine multiple rather than single PPI networks and then benefit from them.

Related work
As mentioned before, we can group the existing approaches to HPO annotations into two categories: feature-based and network-based ones.
Two well-known methods of feature-based approaches are PHENOstruct [9] and Clus-HMC-Ens [18]. Clus-HMC-Ens applies the decision tree ensembles, while PHENOstruct (the extension of GOstruct which was designed to predict GO annotations) relies on the Structural Support Vector Machine (SSVM). Together with HPO annotations (i.e., labels) of each protein, a featurebased method normally accepts feature vectors as the input of a classifier. The trained classifier is then used to make a prediction. The above procedure is the same for both two categories of approaches. Additionally, it is worth noting that PHENOstruct and Clus-HMC-Ens were originally developed for GO but then applied to HPO annotation prediction. In this sense, the difference between HPO annotations and GO annotations has not been fully taken into account by researchers.
Relying on two networks of protein-HPO annotations and the hierarchy of HPO (or Network of HPO, called NHPO) with an optional PPI Network (hereafter PPN), the network-based approaches make predictions. The assumption behind them is that two nodes in a network should share some similarities, particular for those well-connected nodes who have more similarities. In the following, we review the three methods as representatives of network-based approaches, all of which are compared against our proposed approach in the experiments.

Bi-random walk
Bi-Random Walk (BiRW) [19,20] has been demonstrated as a useful method for the bi-network prediction problem. BiRW performs random walks on the Kronecker product graph between PPN and NHPO in a way that they can be combined effectively for the protein-phenotype association prediction. The random walks iteratively performed by BiRW follow the equation: where α > 0 is a decay factor, P and G are the normalized PPN and NHPO matrix, respectively. Y t is the estimation of associations at iteration t, and Y denotes the initial annotations in the training data. By introducing BiRW to capture the circular bigraphs patterns in the networks, the model can unveil phenome-genome associations over time.

Dual label propagation model
The label propagation-based algorithm has been successfully applied to predict phenotype-gene associations in various forms [21,22]. With the following objective function, label propagation assumes that proteins should be assigned to the same label, if they are connected in a PPN: whereS p is a normalized PPN defined asS p = D − 1 2 S p D − 1 2 , and D is a diagonal matrix with the row-sum of S p on the diagonal entries. Equation 2 can be rewritten as follows: where tr(·) denotes the trace of matrix, · F denotes the Frobenius norm, and L S is the normalized graph Laplacian matrix ofS p defined as L S = I −S p .
The Dual Label Propagation model (DLP) [23] extends the label propagation model by adding two smoothness terms. The first term imposes the smoothness in a PPN such that interacting proteins tend to be associated with the same HPO term. The second term imposes the smoothness in NHPO in a way that the connected phenotypes (parent-child pair) are encouraged to be associated with the same protein. The objective function of DLP is given as: where β, γ ≥ 0 are tuning parameters, L S and L G Y encode the PPN and NHPO information, respectively. is the binary indicator matrix that selects only the known associations to be penalized, and denotes Hadamard product (a.k.a entrywise product).

Ontology-guided group lasso
The last method to be reviewed is Ontology-guided Group Lasso (OGL) [24]. It uses an ontology-guided group norm for HPO, rather than the graph regularizer in DLP. By combining label propagation and an ontologyguided group lasso norm derived from the hierarchical structure of HPO, OGL updates estimation, according to the following objective function: where β, γ ≥ 0 are balancing factors. r Y g is the group weight for group g. Y (g)i selects the group members of group g from the i-th column of Y, and the smoothness is imposed through the 2 -norm group lasso ( · 2 ) among the members for the consistent prediction within the group. A notable difference between OGL and our model is that the estimated matrix is not factorized into low-rank matrices.
One of the biggest drawbacks of network-based methods is that data sparseness has a significant impact on the performance. As mentioned before, the current HPO annotations are quite sparse. In addition, all of the network based-methods suffer the heavy computational burden, as they accept a large-scale protein-HPO annotation matrix as an input directly.

Notation
Let Y ∈ {0, 1} N p ×N h be a protein-HPO annotation matrix, where N p and N h are the number of proteins and HPO terms, respectively. If protein i is annotated by an HPO term j, then Y ij = 1, and 0 otherwise. We define S p k (k = 1, 2, · · · , t) be the networks for proteins, namely PPNs, where t is the total number of networks. S p k i,j represents the strength of the relationship between protein i and protein j in the k-th PPN. Similarly, let S h be the network of HPO terms which is generated from an ontology structure and co-annotations, and S h i,j is the similarity value between term i and term j. Our goal is to estimateŶ given Y, S p k and S h .

Our proposed method Preprocessing: generating a network from HPO
The network of HPO terms, or NHPO, is derived by measuring the similarity between two HPO terms in a hierarchy. We adopt the measure proposed in [25]. Having been extensively used in natural language processing, this metric defines the semantic similarity between two labeled nodes by counting the co-occurrence frequency in a corpus.
Specifically for HPO, the semantic similarity between two terms s and t is defined as: where I(s) = log(p(s)) and p(s) = count(s) N p . Here, count(s) denotes the number of proteins annotated by term s and mca(s, t) is given as follows: where A(s, t) represents the set of all common ancestors of s and t.
The weight of the edge between nodes s and t in NHPO is exactly the similarity score. The larger the number of annotated proteins shared by s and t, the higher their similarity score is. It is more likely to happen when the common ancestor of s and t is located closely. This means that S h considers both the co-annotations of two HPO terms and their distance in a hierarchical structure.

Non-negative matrix factorization
The aim of Non-negative Matrix Factorization (NMF) is to find two low-rank matrices with all non-negative elements by approximating the original input matrix. In fact, the latent factors that underlie the interactions are captured. Mathematically, the input matrix Y ∈ R . Then, finding U and V can be done by minimizing the reconstruction error which is defined as: Generally, the 2 (Tikhonov) regularization is imposed to Eq. (7) so as to alleviate overfitting of U and V.
Since there are unknown (missing) entries in Y, we encode the missingness with a masking matrix W ∈ {0, 1} N p ×N h . If the annotation between protein i and HPO term j is missing, we set W ij = 0. Otherwise, we set W ij = 1, meaning that the element Y ij is observed. Accordingly, W is also plugged as an extra input into our model. Together with the 2-norm regularization terms, the objective function is refined as follows: where λ is a regularization coefficient. The unobserved protein-HPO associations are completed by multiplying two factor matrices, or concretely, Y = UV T .

Network regularization
Once we obtain the similarity matrix of HPO, S h , we can regularize V with the help of it. The basic idea is to impose smoothness constraints on the phenotype-side factors; that is where V i is the i-th row vector of V, D h is a diagonal matrix whose diagonals are the node degrees, and L h = D h − S h is the graph Laplacian of S h . Actually, the term is exactly the vanilla graph regularizer. For proteins, multiple PPNs are derived from diverse data sources with heterogeneous properties. In this way, for a collective of PPNs S p k (k = 1, · · · , t), their regularizer is imposed as where L p k = D p k − S p k is the graph Laplacian of S p k , and D p k is the degree matrix. Minimization of graph-based regularization terms will lead to the learned data representations (U and V) that respect the intrinsic geometrical structure of original data spaces (S p k and S h ). Note that such standard graph regularization has already been used in a variety of applications [26].

Model formulation
By combining (8), (9) and (10), our model is formulated as follows: where α and β are regularization coefficients to strike a balance between the reconstruction error and graph smoothness.

Model optimization
Notice that the objective function defined in Eq. (11) is biconvex with respect to U and V. A very regular but effective procedure for fitting is Alternating Least Square (ALS), which alternately optimizes one of the variables by fixing the others as constants until convergence. We first hold U fixed and derive the updating rule of V. The objective function of V can be written as: Accordingly, the derivative of J(V) with respect to V is Taking the Karush-Kuhn-Tucker (KKT) complementary condition, we obtain The multiplicative update rule of V is then: Note that the problem given by (11) is symmetric in terms of U and V. Therefore, the derivation of the updating rule of U is simply the reverse of the above case. Precisely, we have

Training algorithm
We describe the overall framework of HPOAnnotator in Fig. 2. The procedure of our optimization process is presented in Algorithm 1. The optimization was implemented based on the MATLAB code provided by [26].

HPO annotations
Two HPO annotation datasets released by June 2017 and December 2017 were downloaded from the official HPO website (https://hpo.jax.org/). For the sake of brevity, we call them Data-201706 and Data-201712 in the following, respectively. The truepath-rule is applied here to propagate annotations, and only HPO terms with at least one related protein remains. Table 1 lists the statistics of the two datasets.
According to the number of proteins annotated, we separated the HPO terms into five groups: 1 to 10, 11 to 30, 31 to 100, 101 to 300, and more than 300. Figure 3 shows the percentage of HPO terms and corresponding annotations over five groups in Data-201706.

NHPO (Network of HPO)
We downloaded the hierarchical structure of HPO from their official website.

A preliminary test on pairs of two HPO terms in NHPO: the correlation between the number of shared proteins and the average similarity
First, we grouped all pairs of two HPO terms (from NHPO), according to the number of proteins, say M,  shared by the two HPO terms. For each group, we then computed the average similarity score (S h ) by NHPO over those sharing M proteins. Finally, we plotted each group over the two-dimensional space of M× the average similarity score. Figure 4 shows the result. The similarity score is equal to the edge weight of NHPO. This means that this test would be evaluated on the consistency of the similarity with the number of shared proteins from each HPO term pair. There found some correlations between these two, which would be a positive support for using NHPO for HPO annotations.

A preliminary test on pairs of protein-protein edges in a PPN: correlations between the average similarity by a PPN and #shared HPO
Considering the extensiveness, we chosen STRING as the research object. At first step, we grouped all pairs of two proteins, according to the number of their shared HPO terms, denoted as K. For each group, we then computed the average of similarity score (S p ) of STRING PPN over those sharing the same number of HPO terms. Finally, we plotted each group over the two-dimensional space of the average score (similarity) ×K. Figure 5 shows the plotted  Each circle is a pair of two proteins in STRING PPN, with sharing the same numbers of HPO terms, say K. The x-axis is the average similarity score between two proteins over those HPO terms sharing the same K, and the y-axis is K, i.e., the number of shared HPO terms. The red line shows the trend, which is fitted by a polynomial function with the maximum degree of three results. The line in this figure shows that the polynomial trend line is fitted to the distributed points of the twodimensional space. It shows a slightly positive correlation between the number of shared HPO terms and the average similarity score by a PPN. This observation validates the idea that the edges in a PPN may imply that proteins connected by the edges share the same HPO.

Evaluation criteria
The performance is evaluated from three aspects. Annotation-centric measure Each annotation (or a protein-HPO term pair) is viewed as one instance. The models are evaluated using Area Under the receiver operator characteristics Curve (AUC) [31]. Considering the sparseness of protein-HPO association matrix, we measure the Area Under the Precision-Recall curve (AUPR) as well.
Protein-centric measure AUCs (AUPRs) are calculated for each protein based on the corresponding predictive scores by all available HPO terms. Then the computed AUCs (AUPRs) are averaged over all proteins, resulting in micro-AUC (micro-AUPR).
HPO term-centric measure We think that the termcentric measure is important. Typical scientists or biologists focus first on a certain HPO term and are interested in obtaining genes/proteins, which can be annotated by the focused HPO term. The HPO term-centric measure can be computed in a total reverse manner of the proteincentric measure, with the following two steps: 1) AUCs (AUPRs) are first computed for each HPO term; and 2) The computed AUCs (AUPRs) are averaged over all HPO terms, which result in macro-AUC (macro-AUPR).
In addition, we average the computed AUCs (AUPRs) over HPO terms at only leaves of the HPO hierarchical structure. We call the obtained AUC (AUPR) leaf-AUC (leaf-AUPR). We further calculate the macro-AUCs (macro-AUPRs) for each of the five groups, which are generated by focusing on the number of annotations per HPO term (see Fig. 3). In total, (from annotation-, protein-, and HPO term-centric measures) we have the eight criteria to validate the performance.

Experimental procedures Parameter settings
Our approach is compared with three network-based methods: BiRW [20], DLP [23] and OGL [24] as described in related work. Besides, we take Logistic Regression (LR) as a feature-based baseline. Note that LR classifiers are trained on each single HPO term independently, and the features are built by concatenating association scores in PPNs together.
There are several variants of our algorithm by changing the settings of hyper-parameters α and β. We also evaluate each of them as comparison methods. The details are as follows. This setting is in contrast to NMF-PPN. That is, the regularization term of NHPO is kept, while that of PPN is not. For the case of α = 0 and β = 0, there are two another variants depending on whether or not multiple PPNs are utilized.

AiPA: only one PPN is utilized
It is proposed in our previous study [17], which can be regarded as a special case of HPOAnnotator because only single PPN of STRING is exploited.

HPOAnnotator: multiple PPNs are utilized
It is our final model presented in this paper. All four PPNs are used, including STRING, GeneMANIA, BioGRID, and Reactome as described before.

Two evaluation settings
Under two different settings, we validate the performance of the compared methods from two viewpoints:

Cross-validation over Data-201706
We conduct 5×5-fold cross-validation over all annotations on Data-201706. That is, we repeat the following procedure five times: all known annotations are divided randomly into five equal folds. The four folds are for training, while the remaining one is for test. After selecting the test annotation between protein p and HPO term h, all annotations between p and the descendants of term h in the hierarchical structure of HPO are removed from the training data, in order to avoid any overlaps between training data and test data. It means that we predict the annotation of protein p out of all unknown HPO terms, which is a fair and strict evaluation. 2. Independent test by using Data-201712 HPO annotations are incomplete, due to various reasons, such as slow curation. The way of annotations might be changed over time. So we   Table 3 reports the scores of the eight criteria obtained by averaging over 5×5 cross-validation (25 runs in total) on Data-201706. In this experiment, we compare the nine methods in total. In particular, the four are existing methods (LR, BiRW, OGL and DLP), and another five are variants of our model (NMF, NMF-PPN, NMF-NHPO, AiPA and HPOAnnotator). Note that STRING is the only PPN utilized in NMF-PPN. From the table, it clearly shows that our five methods perform better than the four existing methods. For example, our four methods achieve around 0.5 to 0.56 in AUPR, while all the scores by the existing methods are less than 0.1. In fact, our five methods perform better than the existing methods with respect to all of the eight metrics. Thus, their performance differences are very clear. We can conclude that a low-rank approximation is useful for the HPO annotation problem. Furthermore, HPOAnnotator always outperforms other variants in eight conditions among our five methods. This indicates that network information is well incorporated into our formulation. Table 4 lists the AUC scores obtained for five groups divided by the number of annotations. Again, the results reported in these tables demonstrate the same conclusion as that in Table 3. That is, HPOAnnotator outperforms all other methods in all of the cases. A similar trend is also shown in Table 5. In summary, our approach is capable of achieving the best performance for HPO annotations in terms of cross-validation.

Predictive performance in cross-validation on Data-201706
A noteworthy point is that our method works well for the HPO terms with a very small number of annotations, i.e., only one to ten annotations per HPO term. In fact, this situation is usually hard for a low-rank approximation. As HPOAnnotator has achieved the best performance, this implies that a low-rank approximation is useful for all types of groups including HPO terms with a very small number of annotations for HPO annotations.

The effectiveness of individual PPNs in cross-validation on Data-201706
By using NMF-PPN, we perform a set of experiments in order to identify the most effective PPN in terms of HPO predictions. To this end, we perform a series of experiments on NMF-PPN by using a single PPN as its input at a time. NMF-PPN with the four PPNs performs best as reported in Table 6. As shown in Table 6, we can conclude that STRING is the most useful PPN for predicting HPO annotations. By the way, Our model can take advantage of different PPNs to achieve the best performance.

Computation times in cross-validation on Data-201706
The computation (training) times of the eight methods compared in the cross-validation are recorded, where the times are averaged over the total 25 runs (5 ×5 folds). The computation times on the same machine with the same settings are reported in Table 7. From the table, our four models run faster than the compared ones. In fact, they are more than eight times faster than OGL and DLP. The training data is updated periodically, thus the model must be trained by the updated data often. As such, this advantage of our models would make a difference. In addition, OGL and DLP need much more memory spaces than the compared methods.    Table 8 reports AUC obtained by the experiments conducted on independent data for the eight competing methods. Among the three existing methods, DLP achieves the best performance, with AUC of 0.8298. NMF outperforms DLP with AUC of 0.8527, while two variants of NMF with one network regularizer further achieves better performance with AUC of around 0.89. AiPA achieves 0.9187 of AUC with STRING PPN and NHPO. Most importantly, HPOAnnotator archives the best performance, with the AUC of more than 0.92.

Predictive performance in the independent test on Data-201712
As Table 9 reports, seven out of the 30 highest ranked predicted annotations are validated to be true according to Data-201712 which is released later. For example, protein Q02388, encoded by gene COL7A1, is actually annotated by HPO term HP:0001072 (Thickened skin). But we fail to find it in the data released by December 2017. Another example is protein Q9UBX5. According to Data-201706, it has no relationship with HPO term HP:0012638 (Abnormality of nervous system physiology). But this record occurs in the later release of the data.
As the highest-ranked new annotation found by our model, HP:0001072 is known to also annotate another ten proteins, O43897, P07585, P08123, P08253, P12111, P20849, P20908, P25067, P53420, and Q13751, based on Data-201706. We find that their similarity scores with Q02388 in STRING are more than 0.9. It indicates that their interactions between Q02388 and those ten proteins in PPNs imply a high possibility of annotating Q02388 by HP:0001072. In summary, the number of these examples have demonstrated both the effectiveness and necessity of  introducing PPI networks for unknown HPO annotations prediction.

Validating false positives
As mentioned before, seven of the top 30 correct predictions from our model have already been found in the December 2017 release version of HPO annotations. Due to the fact that a curation process on HPO annotations is normally slow, we believe that there may be more false positives among our top ranked predictions. In order to validate our assumption, we first select the rest of the top 10 predictions that have not been found in the December 2017 HPO data. Using a protein name (or its coding gene name) and an HPO term name as a query for online search engines, we then check the relevant literature and diseases for each false prediction. Finally, we manually extract the information from the retrieved papers containing supporting evidence that suggest a particular false positive to be correct in fact. Using this manual process, we find evidence for another two predictions. Table 10 lists the PubMed IDs of the relevant literature, the relevant diseases names, and the detailed evidence for each pair of the found gene/protein-HPO term. The results strongly indicate that the performance of HPOAnnotator is underestimated, which is caused by the incompleteness of the current gold standard.

A typical example of demonstrating the performance of HPOAnnotator
To further demonstrate the performance of our proposed method for predicting HPO annotations, we here present the different predictions made by the four methods for a typical example, protein P23434. As listed in the last row of Table 11, this protein has 10 annotations. It is interesting to note that the number of correctly predicted HPO terms gradually increases from the first row to the fourth row. Again, this indicates that network information is effective for improving the performance of predicting HPO annotations.

Performance comparisons focusing on Organ abnormality
Most of the existing models are evaluated on separate sub-ontologies. However, considering only part of the ontology may lose entire network information. Such information can connect proteins or HPO terms that are even beyond the boundaries of two or more subontologies in the network space. As such, we do not conduct the experiments on separate sub-ontologies. Instead, we focus on the major sub-ontology, Organ abnormality (the part under HP:0000118), with 6370 HPO terms, 3446 proteins and 269420 annotations in total according to Data-201706. A 5×5-fold cross-validation has been conducted by following the same splitting strategy as before. Table 12 reports the scores of the eight evaluation criteria obtained by all compared methods. The results clearly show that the performance differences among the seven cases are subtle. For example, HPOAnnotator achieves the The first three rows of methods with "Organ" are trained by HPO terms on Organ abnormality, while the others with "All" are trained by considering all sub-ontologies. Method performs best in terms of this evaluation metric are in boldface The three rows with "Organ" use only organ abnormality for training, while the others with "All" take all sub-ontologies for training. Method performs best in terms of this evaluation metric are in boldface best performance with respect to all evaluation measure except for leaf-AUC. Comparing NMF-Organ and NMF-PPN-Organ in terms of AUC, we can find that network information can help to improve the performance to a certain extent. Nonetheless, the use of both networks of PPN and NHPO might not be so effective in this scenario. Besides, it seems that the performance improvement is quite limited when we consider the whole ontology rather than individual sub-ontologies. Tables 13 and 14 list the evaluation scores of Macro-AUC and Macro-AUPR over the five HPO term groups, respectively. The trend is similar to that presented in Table 12. Again, the results show no notable difference among the compared methods.

Discussion and Conclusion
In this paper, we have presented an approach that uses a low-rank approximation to solve the problem of the large-scale prediction of HPO annotations for human proteins. In particular, network information is used to regulate such an approximation. The network information can be derived from both sides of annotations, i.e., PPI networks, and a hierarchical structure of an ontology. In essence, we provided a low-rank approximation solution to the optimization problem of matrix factorization with a network-derived regularization. Extensive experiments on the current HPO database have been conducted The three rows with "Organ" use only organ abnormality for training, while the other four rows with "All" take all sub-ontologies for training. Method performs best in terms of this evaluation metric are in boldface to validate the effectiveness of our approach. Experimental results clearly demonstrated the good performance of the proposed method under various settings, including cross-validation, independent test, analysis on the major sub-ontology Organ abnormality, and detailed case studies. The results have validated the good effectiveness as a result of using network information and ontology hierarchical structure as regularization and a low-rank approximation for HPO predictions, even for predictions on HPO terms with a very small number of known annotations. Overall, the four important findings can be concluded from the experimental results: 1) a low-rank approximation works quite well for a large-scale HPO annotations prediction; or more generally, for multi-label classification, even for predicting labels with an extremely small number of labeled instances; 2) a hierarchical ontology structure is very useful as side information for improving the performance of a low-rank approximation; 3) PPI networks from different sources play an important role in predictions; and 4) multiplicative parameter update of a low-rank approximation (matrix factorization) is timeefficient, with around eight times faster than networkbased approaches that need the huge memory because of using the original annotation matrices directly. Authors' contributions JG and SZ jointly contributed to the design of the study. JG designed and implemented the ANMF method, performed the experiments, and drafted the