Figure 1 shows the workflow of ContactPFP. For a query protein sequence, ContactPFP constructs a multiple sequence alignment (MSA) using HHblits (Steinegger et al., 2019), that is, run against the Uniclust30 database (Mirdita et al., 2017) with a parameter set of “-n 3 -id 99 -cov 50 -diff inf”. Then, using the MSA, residue-residue distance prediction (a distance map) is computed for a query using trRosetta (Yang et al., 2020). The predicted distance map of the query is then compared with contact maps of proteins in the reference database using the GR-Align algorithm (Malod-Dognin and Prẑulj, 2014). Since GR-Align compares contact maps of proteins, predicted distance maps were converted to contact maps, or contact graphs, where nodes represent amino acid residues and edges connect residue pairs that are closer than a distance cutoff value. As the distance cutoff values to define a contact, we used 8, 10, and 12 Å between between Cβ atoms. For a given contact distance cutoff (e.g., 8 Å), we define a residue pair as “in contact” if the probability that the pair has the cutoff distance or closer between each other is 0.5 or larger.

The reference database was constructed from Swiss-Prot (The UniProt Consortium, 2021). Sequences shorter than 20 residues and longer than 2000 residues, which made up 1.1% of the all sequences, were excluded. For the remaining 555,378 proteins (98.9%), contact graphs were computed as described above. Using GR-Align, the query contact graph is compared with all contact graphs in the reference database, which are then ranked by graph similarity to the query. In GR-Align, two contact graphs are aligned so that the similarity score of the graphs, which considers graphlet distribution similarity of mapped nodes, is maximized. This graphlet degree similarity can capture the local similarity of contact graphs. A graph similarity score by GR-Align ranges from 0 to 1.0, with 1.0 indicating an exact match in graphlet distributions. Proteins in the reference database that have a similarity score of 0.5 or higher by GR-Align were considered as hits. GO terms from hits were collected and weighted by the sum of the graph similarity scores of hits that have the GO terms. The score of a predicted GO term i is computed as follows: G O S c o r e ( i ) =   ∑ k ∈ P r o t e i n   h i t s   w i t h   G O ( i ) G r a p h _ S i m S c o r e ( k ) (Eq. 1) where k is a hit (i.e., graph similarity score > = 0.5 to the query) that has the GO term i in its annotation. Finally, predicted GO terms for a query is normalized by the highest score among them, so that the most confident GO term has a score of 1.0.

GO terms for proteins in the reference database were compiled from 12 data sources. The primary database used was the UniProtKB/Swiss-Prot. GO terms with IEA (Inferred from Electronic Annotation) evidence code (Boutet et al., 2016) were also included because considering IEA gave a higher function prediction accuracy than excluding them in our previous works (Chitale et al., 2009; Hawkins et al., 2009). In addition to UniProtKB/Swiss-Prot, we integrated annotations from UniPathway (Morgat et al., 2012), TIGRFAMs (Haft et al., 2013), SMART (Letunic et al., 2021), Reactome (Jassal et al., 2020), PROSITE (Sigrist et al., 2013), ProDom (Bru et al., 2005), PRINTS (Attwood et al., 2012), PIRSF (Nikolskaya et al., 2006), Pfam (Finn et al., 2016), InterPro (Finn et al., 2017) and HAMAP (Pedruzzi et al., 2015). This database extension contributed an additional 8,727 GO terms to our reference database. We showed in our previous work that this annotation expansion has a positive effect on function prediction performance (Khan et al., 2015).

We started from representative sequences in UniRef50 (Suzek et al., 2015) (22 October 2019 version). We filter out entries that do not fall within lengths of 100–2000. We only considered the entry names that existed in the 11 November 2019 version of Swiss-Prot. We kept only entries with at least one experimentally verified GO Term in all three categories. To remove the potential redundancy of annotations in the dataset, we only kept one protein from homologous proteins from different organisms. For example, among the two entries, ZRT1_SCHPO and ZRT1_YEAST, which were originally included in the representative sequences, we kept only ZRT1_SCHPO. The ortholog proteins were identified by the common mnemonic protein identification code, e.g.,“ZRTI”. Further, to remove the sequence redundancy, we performed sequence clustering by MMseqs2 (Steinegger and Söding, 2017) with a 25% identity and a 80% coverage (--min-seq-id 0.25 -c 0.8). As a result, we had 9,642 sequences in the benchmark dataset.

To characterize the performance of ContactPFP, we compared it with four sequence-based methods, PSI-BLAST (Altschul et al., 1997), PFP (Hawkins et al., 2009), ESG (Chitale et al., 2009), and Phylo-PFP (Jain and Kihara, 2019). PSI-BLAST is considered the baseline of function prediction methods. We selected PFP, ESG, and Phylo-PFP because they are sequence-based methods that performed well in CAFA challenges (Radivojac et al., 2013; Jiang et al., 2016; Zhou et al., 2019). Our group, who used these three methods, were among the best teams in the series of CAFA challenges. PFP, ESG, and Phylo-PFP use PSI-BLAST search results in different elaborate ways: In PFP, GO terms extracted from PSI-BLAST hits are scored by the sum of the negative logarithm of E-value of the hits. Thus, GO terms from hits with smaller E-values are scored higher. Up to 20,000 hits were considered. In ESG, the top hits of the first PSI-BLAST run are used to perform a second round of database search. In Phylo-PFP, retrieved sequences were ranked by considering both raw E-value and the edge distance on a phylogenetic tree constructed for the sequence. In ESG and Phylo-PFP, raw scores of predicted GO terms are normalized to a range between 0 and 1.0 by the highest score observed for the target protein.

These sequence-based methods identify the query itself as the top hit in a database search. To avoid taking GO terms from the query itself for PFP, ESG, and Phylo-PFP, we removed the query and all hits that had an E-value of 0.0 in the last round of PSI-BLAST before extracting GO terms. These excluded proteins were also removed from the hit list of ContactPFP. For PSI-BLAST, we extracted GO terms from the top 10 hits (except for the query itself and hits with 0 E-value) in the third iteration of a PSI-BLAST run and assigned a score of 1.0 to all the predicted terms.

To start with, we examined how two important hyperparameters in ContactPFP, the definition of residue-residue contacts and choices of top hits from a database search, affect the function prediction accuracy (Figure 2). When constructing a graph from residue distance prediction, a choice of residue distance cutoff needs to be made. A larger distance connects more residue pairs making more edges in a contact graph, while a smaller cutoff would highlight densely connected domains (Yuan et al., 2012). We tested three distance cutoffs between Cβ atoms, 8, 10, and 12 Å, to define residue-residue contacts. The latter parameter, the choice of top hits, decides which proteins from the database search to use for extracting GO terms for annotating the query protein. Increasing this similarity level reduces the number of hits to consider while decreasing the cutoff leads to an increased number of hits.

To examine the effect of the parameters, we performed a four-fold cross validation. In combination with the distance cutoffs, we examined three different ways to select top hits from a search (Figure 2). The Fmax score shown in the panels in Figure 2 are the average values of the four test sets used in the cross validation. In Figure 2A, we used the raw graph similarity score from GR-Align to select hits from a database search. Retrieved proteins in a search that have a similarity score lower than a specified cutoff were discarded. Among the four scores examined, 0.3, 0.5, 0.7, and 0.9, the highest Fmax score of 0.555 was observed when a graph similarity score of 0.7 was used in combination with a distance cutoff of 12 Å. In Figure 2B, we used the top N hits from a database search to extract GO terms regardless of their graph similarity score. The highest Fmax score, 0.638, was achieved when the first two hits were used (i.e., n = 2) with a distance cutoff of 12 Å. In the last panel, Figure 2C, we considered the Z-score of the graph similarity score to select top hits. The highest Fmax score, 0.571, was achieved with a Z-score of 7 using a residue distance cutoff of 12 Å. In each panel in Figure 2, the standard deviations from the four-fold validation were small, 0.002, and the best parameter combinations were consistent across the four-fold.

Overall, the combination of 12 Å for the distance cutoff and using the top 2 hits showed the best performance among the conditions tested. Therefore, we report the results with this condition in the subsequent sections. Regarding the contact distance cutoff, 12 Å showed the best performance in all three panels, which is consistent with what the GR-Align paper reported (Malod-Dognin and Prẑulj, 2014).

We now report the overall performance of ContactPFP in Table 1 in comparison to the other four methods. Two values are reported in Table 1. In terms of the average Fmax score (the left column), ContactPFP’s performance was the second-highest, slightly lower than Phylo-PFP. ESG, PFP, PSI-BLAST followed in this order. Breakdown of the performance in the three GO categories (Table 2) showed essentially the same trend. ContactPFP was the second in Cellular Component, the third in Molecular Function, and the second in the Biological Process.

On the other hand, when predictions given to individual target proteins were compared between two methods (Table 1, the right column), more than half of the proteins (55.6%) had predictions with a higher Fmax score by ContactPFP than Phylo-PFP. ContactPFP also had more wins over ESG, PFP, and PSI-BLAST. Thus, in this head-to-head comparison, ContactPFP was the best. To understand how the performance of the methods differ, we showed the Fmax score of individual target proteins by ContactPFP and each of the other methods in scatter plots (Figure 3). The scores seem not to distribute randomly. Rather, they show an interesting pattern of a “mirror-imaged N-shape”, where there are a substantial number of targets with around 1.0 Fmax as well as other targets with around 0.1 by ContactPFP. This distribution implies that ContactPFP may have both a characteristic strength and weakness when compared with the other methods.

In Table 2, we also presented the Smin score of the five methods. Smin evaluates remaining uncertainty/missing information from predicted GO terms (Jiang et al., 2016). The lower, the better prediction. In terms of Smin, ContactPFP is the best among the all five methods when all GO categories, MF, or BP are considered. ContactPFP was the second in the CC category following Phylo-PFP.

We examined how contact prediction accuracy affects to the GO prediction accuracy in ContactPFP. For this analysis, we used 1,029 targets in the benchmark dataset, which have an experimentally determined protein structure that covers more than 80% of the residues in the target protein. If we consider precision of all predicted contacts (Figure 4A), all the targets fall into contact precision around 0.8 (the average: 0.801), and we found no correlation between the Fmax score. The conclusion was the same when we only considered long-range contacts predicted within the top L/5 scores (Figure 4B); no correlation was observed.

In Figure 4C, we further examined what would happen if we used completely accurate contacts for targets that were taken from experimentally determined structures. Interestingly, the Fmax score was higher when predicted contacts were used. The Fmax score of ContactPFP using predicted contacts and accurate contacts were 0.744 and 0.658, respectively. This is mostly because we use the reference database of proteins with predicted contacts. Similar proteins are likely to have similar predicted contact patterns, either accurate or inaccurate, and the similarity can be captured by contact graph comparison.

We were also curious how ContactPFP performs for proteins that have intrinsic disordered regions (IDRs) because an IDR does not usually form residue contacts. In Figure 4D, we examined Fmax scores of target proteins relative to the fraction of IDRs in a protein. IDRs were predicted with fldpnn (Hu et al., 2021). To make disorder predictions more reliable, we used SPIDER3-single (Heffernan et al., 2018) to predict secondary structure of proteins and only considered residues which were also predicted as loops (class C) SPIDER3-single as the final disorder residues. We did not observe clear correlation between Fmax scores and the fraction of IDRs.

In this section, we discuss cases that illustrate ContactPFP’s performance relative to sequence-based methods.

ContactPFP’s characteristic prediction performance discussed in Figure 3 motivated us to develop ensemble methods. Particularly, the primary focus is to improve predictions for proteins where ContactPFP did not perform well but other conventional methods achieved higher Fmax scores. Combining ContactPFP with other methods also makes sense from a biological point of view because the former uses protein structure information, which is complementary with the latter that uses sequence information.

We constructed ensemble methods of all possible combinations of the methods starting from single methods to a combination of all five methods (Figure 8A; Table 6). When multiple methods are combined, scores of GO terms from the combined methods were simply averaged. The highest average Fmax score, 0.699, was achieved by a combination of three methods, ContactPFP, Phylo-PFP, and PSI-BLAST. Compared with the Fmax of the lone ContactPFP, 0.638, it is a 9.6% improvement. From Figure 4A, we can see that the top methods all include ContactPFP as its ensemble component, which implies it is complementary to the other methods.

Figure 8B shows the Fmax score distribution of the best ensemble method and the five individual methods. ContactPFP’s distribution has a characteristic peak at a low Fmax score of around 0.1. This peak disappeared when combined with other methods, which is the main reason why the performance improved by the ensemble method. In the score comparison of individual proteins (Figure 8C), we can also see that many proteins with a low score around 0.1 by ContactPFP improved by the ensemble method.

One thing which drew our attention is that the top combination included PSI-BLAST, which performed the worst among the five non-ensembled methods. To examine why adding PSI-BLAST improved the performance, we compared the ensemble of ContactPFP with Phylo-PFP with the ensemble of ContactPFP with PSI-BLAST (Figure 8D). We compared PSI-BLAST with Phylo-PFP because the latter was the best single method among the methods we compared (Figure 3A). From this plot, we selected two target proteins where ContactPFP + PSI-BLAST performed significantly better than ContactPFP + Phylo-PFP. APOE_HUMAN, Apolipoprotein E of Human, is one such example. The Fmax score of ContactPFP + Phylo-PFP was 0.234, while that of ContactPFP + PSI-BLAST was 0.801. This protein has 164 GO annotations. By PSI-BLAST, most of the GO terms were found at least once in the top 10 sequences thus all the GO terms had a score of 0.931 (because this is how GO terms are scored in PSI-BLAST). These GO terms were found also in the sequences retrieved by Phylo-PFP and ContactPFP. However, the GO terms were found infrequently in the retrieved sequences and thus received a low score. The same mechanism was observed for DNJA3_HUMAN, DnaJ homolog subfamily A member 3, which had an Fmax score of 0.213 by ContactPFP + Phylo-PFP and 0.988 by ContactPFP + PSI-BLAST.

In Figure 9 we show the computation time of ContactPFP. In our computational environment, the entire ContactPFP pipeline took approximately 15 min for a 500 residue-long protein. Figure 8 also shows the breakdown of the time needed for five steps in ContactPFP. The computational time for running trRosetta and reference database search by GR-Align grows as the protein length increases. Particularly, the computational time for using trRosetta grows sharply, and it exceeds the time for the reference database search when the query protein is longer than 500 residues.