Predicting phenotypic traits of prokaryotes from protein domain frequencies

Investigation of discriminative domain families

In order to assess the biological relevance of the discriminative protein domain features, we investigate the 50 most discriminative phenotype-specific domain families within a biological context. For brevity, we restrict the analysis to the two phenotype categories "Endospores" (spore formation) and "Motility". However, the investigation can be performed in a similar fashion for the "Gram stain" and "Oxygen Requirement" phenotypes using the ranked lists of discriminative domain families provided in additional file 2. The features are based on the discriminative model that has been learned from all training organisms using the parameter associated with the best validation performance. We also discuss the results from the phylogenetic clustering of these families, which may be used to identify functional modules of domains (for details see section "Methods"). The clustering dendrograms associated with all phenotypes can be found in additional file 3.

Endospores

Endospores are dormant resistance forms produced by some gram-positive bacteria by the formation of a thick internal wall that encloses the DNA and a part of the cytoplasm. For reasons of brevity, table 4 only shows the 10 most discriminative domain families associated with the phenotype category "Endospores", a list with the 50 most discriminative domains can be found in additional file 2. Of the top 10 (50) domains, 4 (12) have the terms "spore" or "sporulation" in their description line. Additional domains are directly involved in sporulation, e.g. the Cell wall hydrolase, the Coat F domain, the Germination protease, or the YabG peptidase U57 domain, which is present in a sporulation-specific protease of several spore coat proteins [18]. Other domains such as YabP have been associated with the sporulation process through mutagenesis analyses. Thus, the discriminative domain set for prokaryotes associated with the term "Endospore formation" is enriched in domains involved in the sporulation process.

rank	weight	# groups	Pfam-ID	Pfam description
1	0.008	3	PF03419	Sporulation factor SpoIIGA
2	0.007	5	PF07486	Cell Wall Hydrolase
3	0.007	1	PF06686	Stage III sporulation protein AC (SpoIIIAC)
4	0.007	2	PF00269	Small, acid-soluble spore proteins, alpha/beta type
5	0.007	1	PF07873	YabP family
6	0.007	1	PF09555	Stage III sporulation protein AD (spore_III_AD)
7	0.007	3	PF00407	Pathogenesis-related protein Bet v I family
8	0.007	4	PF00876	Innexin
9	0.007	6	PF04672	Protein of unknown function (DUF574)
10	0.006	1	PF04647	Accessory gene regulator B

List of the 10 most (positive) discriminative domain families associated with the RLSC model for the phenotype category "Endospores" (see also section "Methods"). The first column indicates the rank, the second column shows the discriminative model weight. The third column denotes the phylogenetic width of a particular domain family, i.e. the number of taxonomic groups at phylum level in which the family occurs. The fourth and fifth column correspond to the Pfam ID and family description associated with a particular domain family. The table with the 50 most positively and negatively discriminative domain families can be found in additional file 2.

Figure 2 shows the dendrogram resulting from phylogenetic clustering of the 50 most discriminative domain families. In the upper part, a cluster of Pfam domains associated with sporulation can easily be identified (colored green). Interestingly, also three domains of unknown function ("DUFs") are localized in that cluster (DUF1659, DUF1540, DUF1429). This indicates that the proteins associated with these families may directly be involved in the sporulation process or act together with sporulation proteins in the same pathway. In fact, the domain family DUF1429 (PF07241) located at the lower part of the cluster is considered a "dead" family in the current version (24.0) of the Pfam database and has been merged into the family YabP (PF07873), which is associated with sporulation.

Detecting annotation errors

A possible use of phenotype prediction approaches is the completion and correction of existing annotation. In order to identify erroneous NCBI annotation associated with the organisms, we compared our prediction results with the phenotype annotation and verified the candidates via literature search. We found several cases of possible annotation errors for different categories by verifying putatively false predictions on the test set that do not agree with the existing annotation. For instance, for the "Motility" category we identified two organisms (Burkholderia glumae BGR1, Bacillus cereus 03BB102) with negative annotation but evidence for being motile [20, 21]. Similarly, for the category "Gram stain" we found another two possible erroneous annotations (Eubacterium eligens ATCC 27750, Geodermatophilus obscurus DSM 43160) with phenotype annotation in the literature differing from that of the NCBI [22, 23]. Finally, we identified two organisms that most likely have been assigned to the wrong respiratory subcategory (Dethiosulfovibrio peptidovorans DSM 11002, Acidimicrobium ferrooxidans DSM 10331[24, 25]). In order to encourage the systematic correction of erroneous existing annotation, we provide phenotype-specific lists with prediction results from our approach for organisms in the independent test in additional file 4.

Robustness and speed

Our validation (test) data set consists of 795 (106) finished and 237 (337) unfinished genomes, respectively. The total number of examples used for model training and validation is substantially higher as compared to related work of e.g. Slonim et al. [6] (59 genomes) or Kastenmüller et al. [8] (266 genomes). We think that the robustness of models, which is indicated by the good generalization performance, can partly be attributed to the unprecedented size of the data set. Furthermore, the high prediction accurracy on the test set, which contains many yet unfinished genomes, suggests that even in early stages of genome projects reliable statements about an organism's phenotype can be derived from the predicted domain profile.

In contrast to the approach presented in [8], our method does not involve an intermediate mapping of proteins to pathways. Instead, the phenotype prediction is directly performed on the basis of estimated Pfam domain profiles. The computation time that is necessary for inference of the domain profiles depends on the domain detection method and the size of an organism's proteome and ranges from a few seconds to several hours [15]. The actual phenotype prediction can be efficiently performed by means of the dot product of the precomputed domain profile and the discriminative weight vector in domain profile space. This only requires storage of the 10797-dimensional discriminant and ≈ 0.05 milliseconds for calculation of the prediction score.

Dataset

We downloaded all prokaryotic proteomes from the integr8 FTP site ftp://ftp.ebi.ac.uk/pub/databases/integr8/fasta/proteomes/ on May 12, 2009 and selected all NCBI annotated organisms which could be uniquely associated with an integr8 proteome file via the integr8 proteome report ftp://ftp.ebi.ac.uk/pub/databases/integr8/proteome_report.txt. The annotation of these 1032 organisms for the phenotypes "Endospores", "Gram stain", "Motility" and "Oxygen Requirement" was extracted from the respective columns of the table on the NCBI prokaryotic genome project web site ftp://ftp.ncbi.nlm.nih.gov/genomes/genomeprj/lproks_0.txt. While these proteomes were used for the validation procedure (see below), we downloaded additional 443 proteomes and the associated annotation from the abovementioned websites on October 10, 2009. The complete list of organisms and their phenotype annotation can be found in additional file 1.

Construction of Pfam domain and COG profiles

In this work, the Pfam domain profile associated with a particular organism is a d = 10797 dimensional domain occurrence vector according to the maximum family index of PF10797 in version 23.0, section "A" of the Pfam database [13]. This vector has nonzero entries only for dimensions associated with the Pfam domains that occur in the organism's genome. To calculate the occurrences, we applied a Pfam domain detection to each protein in each genome using the UFO method [15]. The resulting vector x of absolute domain occurrence counts is normalized to relative domain frequencies such that $x_f\in {[0,1]}^d,{\displaystyle {\sum }_{i=1}^d{x}_i}=1$

Furthermore, each dimension of the domain feature space with a standard deviation different from zero is normalized to unit standard deviation. The data matrices with UFO counts for all organisms of the validation and test set along with the organism names and relevant NCBI phenotype annotation in comma separated value (CSV) format can be found in additional file 5.

For comparison, we also calculated organism-specific profiles of full gene frequencies in terms of clusters of orthologous genes (COGs, [17]) detected in the organisms' protein sequences. For COG detection, we downloaded the COG database in RPS-BLAST (Reverse PSI-BLAST, [26]) format from ftp://ftp.ncbi.nih.gov/pub/mmdb/cdd/little_endian/Cog_LE.tar.gz and evaluated the organisms' proteomes against this database using the default e-value threshold of 0.01. The 5665 COG clusters give rise to 5665-dimensional organism-specific feature vectors, which we normalized as described above.

Phenotype prediction on domain profiles

To predict microbial phenotypes from organism-specific domain profiles, we use a supervised classification approach. For this purpose, we divide the task of predicting four phenotype categories into four binary classification problems according to the discrimination of different phenotype realizations. Here, organisms that are annotated as "yes" ("Endospores"), "+" ("Gram stain"), "motile" ("Motility") or "aerobic" ("Oxygen Requirement") are considered as positive examples and organisms that are annotated "no"/"-"/"non-motile"/"anaerobic" are used as negative examples of the respective two-class problem. Domain profiles of organisms that do not have any annotation for a particular phenotype realization are not considered for evaluation. Furthermore, organisms with an identical UFO Pfam domain profile were reduced to a single representative organism. Thus, the total number of examples as well as the number of positive and negative examples vary across the different phenotype categories. An overview of the number of examples associated with each phenotype and the phylogenetic distributions of all phenotypes in terms of histograms can be found in additional file 6.

where B is a matrix that contains the inverse elements of the vector b as diagonal elements. In case of multi-class learning problems with M classes, the label vector y has to be replaced by a matrix Y = [z₁,..., z _M ], z _{i, j} ∈ {0, 1}, where z _i is a class-specific vector with non-zero values for examples belonging to class i. The minimization then yields a matrix W = [w₁,..., w _M ] of M class-specific discriminative weight vectors.

To evaluate the influence of the regularization parameter λ, we randomly divided the data set into 20 partitions with 70% training and 30% validation examples, respectively. Using these partitions, we computed the average area under curve with respect to the precision-recall characteristics (PRC, [16]) over all partitions to determine the best parameter λ = {10 ^m |m = -5, -4,..., 5}. Finally, we tested the prediction performance of our method in terms of the harmonic mean (also known as F₁-measure) $harm=\frac{2*sens*spec}{sens+spec}$, which combines sensitivity $(sens=\frac{TP}{TP+FN})$ and specificity $(spec=\frac{TP}{TP+FP})$. To estimate the generalization performance of our approach on the set of 443 test genomes, we evaluated the discriminative model associated with the highest validation performance.

The kernel-based RLSC model associated with a phenotype classification problem is represented by a vector of N organism-specific weights. For fast prediction of phenotypes, the discriminant in the original feature space of Pfam protein domain profiles can be calculated by a linear combination of the learned organism-specific weights and the domain profiles in X. The phenotype prediction for a newly sequenced organism then only requires the construction of the organism's Pfam domain profile and the computation of the dot product of the feature space discriminant and the domain profile.

The feature space discriminant also allows to inspect the learned discriminative features in terms of phenotype-specific Pfam domain families. For each phenotype prediction problem, we assemble two ranked lists of domain families that are associated with the 50 largest positive and negative weights in the profile space discriminant, respectively. This allows to identify indicative and counter-indicative domain families directly.

For further interpretation of the abovementioned lists and for identification of functional modules, we also calculate clustering dendrograms for the 50 most discriminative domain families based on their phylogenetic profile, i.e. the absence/presence pattern of a particular domain across all organisms. To calculate the dendrograms, we apply the Matlab^® function 'dendrogram' to the domain phylogenetic profiles using average linkage (UPGMA) and the profile correlation as a distance measure. Dendrogram branch colors are calculated by the 'dendrogram' function using the default color threshold at 70% linkage branch length.

Comparison with pathway-based prediction

For comparison to a method that can in principle be used for microbial phenotype prediction, we evaluate our approach using the dataset presented in [8]. Because of the different set of phenotypic traits used in the two works, only the two phenotypes "Gram stain" and "Oxygen Requirement" can be directly compared. As in [8], we selected the organisms associated with the phenotypical traits "aerobic"/"anaerobic" and "gram positive"/"gram negative" as the respective positive and negative examples of the phenotype category. In Kastenmüller et al. [8] the commercial Pedant database [30] has been used for collection of genome features and phenotype annotation. As a consequence, not all organisms present in the Pedant database can be found in the NCBI prokaryote genome project. For our evaluation, we used the data of all organisms whose taxonomy ID could be mapped to a NCBI taxonomy ID. In total, we used 166/92 out of 201/113 possible organisms for the phenotype categories "Gram stain"/"Oxygen Requirement", respectively. However, our inspection of the list of all organisms revealed that most of the missing organisms have identical species names when compared with the present organisms, thus the performance values we obtain here can be expected to provide lower bounds.

For direct comparison of the results, we measure the prediction performance in terms of the area under ROC curve [16] and the product of sensitivity and specificity as in [8]. In the original work, a ten-fold cross validation procedure has been performed. To get a reliable estimate of the prediction performance, we repeated the cross validation 100 times using random partitions. Furthermore, we compared our approach to the best of the four classification methods used in [8] as indicated by the classification quality diagrams in additional file 2 of the original work.