An en masse phenotype and function prediction system for Mus musculus

Compilation of data

Data of a large number of distinct types were collected and organized for both phenotype and GO term prediction. Protein domain annotations, protein-protein interactions, expression data, disease data, and phylogenetic profile data for 21,603 M. musculus genes were taken from the dataset established for the MouseFunc project [14]. Expanded and more recent (relative to the MouseFunc data [14]) sets of mammalian phenotype and GO term annotations for the same mouse genes were acquired from Mouse Genome Informatics (MGI) [17]. Detailed descriptions of the data organization and pre-processing are given in the Materials and methods (see below).

A combined learning approach for annotation prediction

A machine-learning approach was used here for function and phenotype prediction. We used a weighted combination of two distinct models, rather than adopting a single family of models for prediction. For clarity, a brief description of the general method is described here, with further details available in this issue [15].

First, we applied a gene-centric 'guilt-by-profiling' approach. Here, properties associated with single genes (for example, matches to specific protein domain patterns) were used to infer various annotations for those genes. In this case, a distinct random forest classifier [18] was constructed for each GO term or phenotype. Second, a paired-gene approach was applied ('guilt-by-association'). Pairs of genes are connected by weighted edges, with weights representing belief that the two genes are functionally linked (that is, a functional-linkage graph [19]). Genes with strong functional links tend to share annotations, and in cases where one member of a functionally linked pair has an annotation and the other does not, a new annotation can be inferred.

Annotation predictions made by each of these two component classifiers were then combined using a regression model that maximizes a chosen performance measure, exploiting the best features of both predictive classifiers. A brief description is given in the Materials and methods (see below) and in more detail in an accompanying paper describing the methodology and its application to S. cerevisiae gene function [15].

Genome-wide application of a combined learner to M. musculus

Predictive models for 2,938 GO and 3,914 mammalian phenotype terms were made; 21,603 M. musculus genes were given quantitative scores for each possible annotation of each term (for a total of 21,603 × [2,938 + 3,914] ≈ 148 × 10⁶ prediction scores computed).

Quantitative predictive performance estimates

To evaluate the participating classifiers, the MouseFunc project employed the area under the receiver operator characteristic (ROC) curve (AUC-ROC) and precision at specified levels of recall (precision is the fraction of predictions [scores above a given threshold] that are true, while recall is the fraction of true annotations recovered above some score threshold). We use the same general approach here, but include an additional measure: 'mean average precision' (MAP). An ROC curve indicates the relationship between true positives and false positives as the score threshold for calling a prediction is varied. Also used are points along the precision-recall curve (as the threshold varies), with MAP being the mean precision obtained at all distinct recall levels. In the case of a tie among scores, the average precision over all permutations of the response variable (the GO or phenotype term annotations) is taken as the precision for that level of recall. The MAP has been identified as a good alternative to other precision-recall based statistics. Most notably, it has been deemed more useful than area under the precision-recall curve (AUC-PR) as a measure of comparison between classifiers [20]. We also include aggregate performance measures for categories such that scores for all terms within a single category are pooled.

Table 2 gives - for each category - the aggregate precision for the combined classifier across varying levels of recall, as well as the AUC-ROC. Precision at low recall is usually the driving criterion for biologists, who generally want to have high precision (accuracy) when examining a 'short list' of candidate genes for association with a particular function or phenotype.

	P01R*	P05R	P10R	P20R	P50R	AUC-ROC
BP [3,10]	0.41	0.23	0.17	0.08	0.002	0.75
BP [11,30]	0.41	0.28	0.22	0.14	0.01	0.79
BP [31,100]	0.57	0.51	0.42	0.30	0.06	0.75
BP [101,300]	0.60	0.49	0.43	0.32	0.09	0.85
CC [3,10]	0.78	0.43	0.45	0.24	0.02	0.85
CC [11,30]	0.72	0.54	0.49	0.38	0.04	0.86
CC [31,100]	0.71	0.65	0.51	0.36	0.04	0.86
CC [101,300]	0.92	0.72	0.59	0.46	0.12	0.87
MF [3,10]	0.65	0.52	0.52	0.48	0.25	0.88
MF [11,30]	0.70	0.68	0.62	0.56	0.34	0.90
MF [31,100]	0.84	0.74	0.66	0.60	0.43	0.93
MF [101,300]	0.86	0.70	0.69	0.61	0.47	0.92
Pheno [3,10]	0.10	0.07	0.03	0.003	0.0006	0.70
Pheno [11,30]	0.20	0.13	0.07	0.02	0.003	0.75
Pheno [31,100]	0.31	0.20	0.14	0.07	0.009	0.78
Pheno [101,300]	0.44	0.33	0.25	0.17	0.04	0.78

*P01R is the precision at 1% recall, and the other headers follow this pattern. AUC, area under the curve; BP, biological process; CC, cellular component; MF, molecular function; Pheno, phenotype; ROC, receiver operating characteristic.

For GO term prediction, mean AUC-ROC for the combined classifier exceeds 0.8 for 9 out of 12 categories (with 0.5 AUC-ROC being expected for random predictions). At 1% recall, the precision ranges from 41% to 92% across the 12 categories, with MF and CC terms being easier to predict for than BP terms, and with difficulty in prediction increasing as the existing annotation count for these terms decreases (due to limited training data availability and a lower prior probability of success due to low prevalence of existing annotations). The performance levels for each category as the threshold varies can be seen in Figures 1 and 2.

Phenotype annotations prove systematically more difficult to classify than GO term annotations (see Discussion, below). The AUC-ROC for each pooled phenotype category exceeds 0.7, with precision at 1% recall ranging between 10% and 44% (again with the observation that terms with lower annotation count have lower predictive accuracy). Performance characteristics along the threshold range are depicted in Figures 3 and 4.

Examination of top novel predictions

Cross-validation and out-of-bag [18] performance estimates are useful to assess the ability of the classifier to recover the 'gold positives' not used in the training process and, thus, the generalizable performance of the classifier. However, the true potential of such classifiers to impact biology is in their ability to make novel annotations, even for genes used in the training process. Thus, high-scoring 'incorrect' predictions are of extreme interest because they provide novel hypotheses about gene function and may serve to guide choice of experiments. (We note that some predictions may only be 'incorrect' in the sense that cross-validation or held-out test set performance measures consider such cases as 'false positives' because the gene has not yet been annotated with the property, even though such a prediction may in fact be correct and thus 'novel'.) The MouseFunc project introduced a method for prospectively assessing the ability of classifiers to make novel predictions [14]. Training data for MouseFunc reflected the state of annotation knowledge in February 2006. When scores were submitted in October 2006, the classifiers' ability to recover held-out test data was assessed, but in addition the classifiers' ability to predict for annotations that were newly assigned between February and October was measured. The methods in the MouseFunc project (group 'G') are fundamentally the same as those applied here; thus, we expect approximately the same level of novel term prediction performance while using a newer set of GO data for training. It should be noted that the performance reported in that work can be viewed as a pessimistic measure of the actual performance, as only a matter of months elapsed before assessing novel predictions, and many of the apparently false-positive predictions at the end of that period may ultimately be proven correct by future experiments.

To gain intuition from specific examples, we examined some of the most interesting novel predictions within the literature. A gene/term prediction is novel and was considered interesting if it was not currently annotated in the reference database and if there did not exist any current annotation involving the gene and any non-root ancestor of the term. Otherwise, we consider such predictions to be 'refinements' of existing annotations. Within each category, the top interesting gene/term combinations are scanned in decreasing order of confidence (that is, prediction score) and - to avoid over-weighting particular genes or terms - a further filtering step is employed to limit each gene and term to a single appearance in the list to be followed up in depth. This filter is described in detail in the Materials and methods (see below).

Availability of predictions

All predictions have been made available through the worldwide web through a simple database gateway [29].

Success of a combination-of-learners approach

'Guilt-by-profiling' and 'guilt-by-association' methods have both been previously employed for functional prediction. The combination of these approaches has been applied to the unicellular eukaryote S. cerevisiae [13, 15]. Here we show that an approach that combines classifiers of both types by logistic regression is superior to that of either base classifier alone in the mammal M. musculus for both function and phenotype. A detailed discussion of this combination approach is provided in [15]. Briefly, a single free parameter α weighting the contribution of each base classifier was chosen to maximize AUC-PR. To avoid over-fitting, a single α was found for each pooled category. Performance evaluation of the algorithm as applied to M. musculus function and phenotype is described below.

As observed in [15], each 'base' classifier type (for example, functional linkage decision trees) has unique strengths in different categories and across different performance measures. The functional linkage approach ('guilt-by-association') tends to be the favored base classifier when evaluated by ROC-based measures (Figures 5 and 6), but the 'guilt-by-profiling' method excels when precision is the performance measure of choice (seen in Figure 7 and to a lesser extent in Figure 8). We believe that for many researchers the precision at low recall levels is of greatest importance - and is ultimately how the optimized combination method was guided in this study, aiding in the prioritization of a few high-scoring results for a laborious literature-based follow up.

As seen with S. cerevisiae [15], different base classifiers may each excel for some categories, or within a bounded specificity range of a single category (for example, seen in some cases in Figures 1 and 3). For prediction of M. musculus GO terms, Figure 2 shows that CC terms seem to benefit the most from using the functional linkage approach and that MF terms do not gain much from the integration of the functional linkage data. In contrast, MF annotation predictions (often relying on specific domain features - inherently gene-centric data - as primary predictors) tend to rely on the output of the random forests approach. By tailoring the data most appropriate for each style of base classifier, the combination approach leads to greater overall performance [15].

Figures 1, 3, 7, and 8 show the increase in precision across the entire recall range (and as a single MAP score) when using the combined scores relative to performance of each base classifier alone. Because we optimize the combination such that AUC-PR is maximized, we also checked to see whether the MAP measure for individual terms also increased. Figures 7 and 8 show that the mean of the MAP within each category is increased using the combined approach. This is consistent with results for S. cerevisiae function prediction [15] but had not been demonstrated for mammals, for phenotype prediction, or using the MAP measure. The combined classifier's performance in a precision-recall statistic does not guarantee improvement in other performance measures. However, the combination optimized for AUC-PR also led to AUC-ROC improvements across most GO categories (Figure 2) and had little impact for phenotype categories (Figure 6). We believe that predictions with the most utility are those with high precision, which tend to be at low recall. Thus, we conclude that the substantial benefits gained in precision using our regression-based precision-optimized combination easily justify any potential drop in AUC-ROC measures. The ability to increase predictive power using a combination of two classifiers leads to the natural question of, 'why stop at 2?' Indeed, one can combine in a similar fashion many different classifiers, with the number of degrees of freedom in the optimization problem increasing linearly with the number of base classifiers. Comparison of the optimized multiplicative coefficients might then be used to rate the relative contribution of each classifier to the final combined model, while simultaneously providing the best performance possible under the constraints of the method for combination chosen.

A method for characterizing high-scoring predictions

As discussed in [15], random forests offer an advantage over some alternative machine-learning methods (for example, support vector machines, or monolithic classification trees) in that they provide a reliable variable importance measure. This provides valuable insight into why a prediction was made, offering the researcher a first step in understanding the biochemical principles underlying the annotation. As an example, a class (i) prediction associating the gene Plscr1 (phospholipid scramblase 1) and the GO term 'cholesterol homeostasis' used the Interpro-catalogued protein domain pattern 'Proteinase inhibitor, hirudin/antistatin' as the primary predictor (statins play a role in cholesterol regulation through inhibition of HMG-CoA reductase activity).

Data sources

For both GO and mammalian phenotype term annotation prediction, our training data consisted primarily of the datasets assembled and integrated for the MouseFunc project. While much of the data originates from outside of the mouse community, the curation and data integration (for example, mapping between human interactions and their mouse equivalents) was primarily taken from the Mouse Genome Database [4] as part of the organization of the MouseFunc project. Cases where data used here differed from that used in the MouseFunc project are made clear below.

Training data organization for functional linkage approach

Training data for the machine learning algorithms were organized generally as described in [15], with descriptions specific to M. musculus data manipulation given below.

The Pearson product-moment correlation coefficient of expression data was computed for each gene-pair (for each of the three data sets). These continuous coefficients were then binned into five groupings, E_0.5, E_0.6, E_0.7, E_0.8, and E_0.9; each group defined as:

E _x ={(a, b)|ρ₁(a, b) ≥ x}

where a and b are genes, and ρ₁ is the correlation coefficient function. This resulted in 15 binary attributes describing each gene-pair. The protein-protein interaction data were already provided in gene-pair format, with each positive edge resulting in a positive gene-pair for a single binary matrix column representation.

Those gene-pairs (a, b) having ρ₂(a, b) ≥ 0.9 were assigned membership in a single set. Each domain dataset (Pfam and Interpro) was processed using this method, resulting in two binary variables.

The phenotype and OMIM disease data were given the same treatment as that for the protein domain data, resulting in two additional binary variables. The phylogenetic data sets were processed identically to the protein domain data, except that only three final Jaccard thresholds were used (0.7, 0.8, and 0.9), resulting in 6 additional binary variables (3 each for Biomart and Inparanoid data).

Thus, a total of 26 distinct binary variables describing each gene-pair were used by the functional linkage classifier. For phenotype prediction, eight additional variables were included representing co-annotation of terms within the MF and CC annotation sets (four groupings each based on annotation counts in ranges [3,10], [11,30], [31,100], or [101,300]), leading to 34 binary variables. We excluded terms within the BP branch to avoid potential circularity because GO BP terms are often tightly related to phenotypes. The data were then organized into a matrix with genes as rows and the 34 predictive variables as columns.

Training data organization for random forests approach

For the random forest base classifiers, binary gene-centric data were represented in a binary matrix form (each row representing one gene, each column representing a property associated with that gene). When predicting GO terms, the training matrix was composed of the protein domain data, the phenotype data (top level only), the phylogenetic profile data, and the disease data (for a total of 5,405 + 3,133 + 33 + 18 + 24 + 2,488 = 11,101 columns). For phenotype term prediction, the training matrix consisted of the same data used for GO prediction as above except without the phenotype data (for a total of 11,068 columns).

Probabilistic decision trees for functional linkage approach

To generate functional linkage graphs, a single probabilistic decision tree [44] was built for each category (12 for GO prediction, 4 for phenotype prediction). For each tree, the response variable was whether or not the gene-pair shared an annotation of a specific term within that category or any category in the same branch but with smaller annotation count. This response variable is thus used as a proxy for the general concept of two genes being 'functionally linked'. For GO term prediction, training examples (gene-pairs) were filtered such that only those genes with some existing GO annotation were included in the gene set, leading to $\left(\begin{array}{c}10542\\ 2\end{array}\right)$ unique gene-pairs. No such filter was performed for the four phenotype functional linkage trees (giving $\left(\begin{array}{c}21603\\ 2\end{array}\right)$ gene-pairs). Decision tree training and classification was then performed as described in [15].

Random forests approach

Random forests are used as described in [15], with parameter choices and organism-specific details as listed below. For each GO term prediction a forest of 300 trees was built, with a column sample size chosen as v/20 where v is the number of remaining variables to be sampled at a given node. For the phenotype predictions, 200 trees/forest and the same v/20 were used as parameters. Each term (2,938 GO and 3,914 phenotype) corresponds to a separate random forest model.

Logistic combination of methods

Calibration of scores

where count(t) is the number of annotations for term t (that is, the average scaled score will equal the prior probability of positive annotation for term t).