Measuring pathway database coverage of the phosphoproteome

Systematic literature review of public databases

There is an abundance of pathway databases currently available for public use in interpreting biological data (Kanehisa & Goto, 2000; Joshi-Tope et al., 2005; Slenter et al., 2018; Zhou et al., 2012; Wishart et al., 2006; Pratt et al., 2015; Szklarczyk et al., 2019; Cowley et al., 2012). However, they tend to differ from each other in several properties. To address this, we systematically compared 29 databases and produced a selection criterion for inclusion in our analysis (Table S1, Supporting Information).

To be included in this comparison, each database had to fit the following criteria:

• Have phosphorylation information

• Have a connected signalling network component

• Be updated yearly, and if not must have been widely used and be widely cited (1,000 minimum)

• Freely accessible

• Downloadable

The pathway databases which met these criteria were HPRD, BioGRID, SIGNOR, Reactome, KEGG, PhosphoSitePlus and WikiPathways.

Data retrieval and processing

Since many of these databases were represented in differing formats and some with unique identifiers, several processing steps were conducted as needed per database to extract protein UniProt accessions and phosphosites as outlined in Fig. S1, Supporting Information. Dates of database access are outlined in Table S2, Supporting Information. For those databases available in a BioPax format, they were first converted into a Neo4j graph database object (Neo4j, 2020) (Schema described in Fig. S2, Supporting Information). The Neo4j graph database implementations were built and used to allow for quick and easy extraction of relevant information as described in more detail below.

Protein

For IMEX, qPhos and, UniProt and both Gene Ontology (GO) categories the data was available in TSV format with UniProt accessions easily extractable. Each database downloaded in the BioPax format (Reactome and PhosphoSitePlus) was converted into a Neo4j graph database in Java, where the list of UniProt accessions were extracted. SIGNOR’s API was used to download the list of UniProt IDs involved in the signalling network. BioGRID’s PPI file was downloaded in TSV format directly from the website with RefSeq IDs from the signalling portion of the database being taken. HPRD legacy files were downloaded from the Integrated Network and Dynamical Reasoning Assembler (INDRA) python package documentation (Gyori et al., 2017). For HPRD, RefSeq identifiers were taken, once again only from the signalling portion of the database. We consider the interaction network portion of HPRD and BioGRID, as a source of connected proteins with potential signalling function. For both KEGG and WikiPathways only the identifiers belonging to a pathway were downloaded as we were only interested in proteins involved in signalling. The KEGG database contained only KEGG identifiers, and the WikiPathways download contained only Entrez IDs. All sets of identifiers were converted to UniProt accessions. No identifier type mapped cleanly to UniProt IDs, so in cases where more than one UniProt accession was found, only the reviewed UniProt accession was taken. If there were multiple reviewed UniProt accessions found, the UniProt accession with the longest protein sequence was selected. If no reviewed UniProt accession was found, the unreviewed UniProt accession with the longest sequence was taken. This was done to maintain a one-to-one ratio of alternate identifiers to UniProt accessions. Because some resources annotate pathways with gene level identifiers, we remove annotated isoforms from other resources in order to maintain a fair comparison.

Once the list of UniProt accessions were obtained for each resource, each was updated to the latest version of UniProt (Jan. 2021). In this instance, ‘update’ refers to converting all secondary accessions to primary accessions, removing all non-human identifiers and removing identifiers that are classified as obsolete in UniProt.

Phosphosites

In Reactome and PhosphoSitePlus, the phosphosites were extracted from the neo4j graph databases generated for extracting the UniProt accessions mentioned above. For qPhos, the set of all phosphosites were extracted from the main data file downloaded. Phosphosites from HPRD were extracted from the separate PTM file and connected to UniProt IDs converted from RefSeq IDs with only phosphosites belonging to the signalling portion of the database taken. BioGRID’s phosphosites were also extracted from a separate PTM file and connected to UniProt IDs converted from RefSeq IDs only from the signalling portion of the database. For SIGNOR, the human PTM file was downloaded directly from the website and processed. Phosphosites from UniProt were extracted directly from the January 2021 downloaded TSV file. Only phosphosites from reviewed UniProt accessions were used in this study. Upon further investigation, phosphorylations listed in KEGG had no specific site annotations. WikiPathways contained multiple phosphosite annotations in multiple formats, making it impossible to extract them all. As with above, since each phosphosite was associated with a single UniProt accession, the UniProt accessions were updated to reflect only human and current UniProt accessions and isoforms were removed.

Consistency analysis

To perform the consistency analysis the neo4j graph database version of Reactome and PhosphoSitePlus were used. For SIGNOR, HPRD and BioGRID the list of unique phosphosites retrieved from each respective database file was used. In each database, all PTM’s belonging to a unique UniProt id were gathered into a set. A phosphosite in this manuscript is defined as a single amino acid reported to be phosphorylated and its position as reported in the database. We do not take into surrounding sequence or attempt to standardise or harmonize this information across databases beyond what is described in the methods. Then, for each UniProt ID the canonical sequence and all isoform sequences were retrieved from UniProt. Next, for each PTM annotation, if the amino acid at the annotated position was the same in the canonical UniProt sequence it was called a match; for example, if Reactome contains an annotation on Protein X that indicates a phosphorylation on a threonine at location 128 and UniProt has a threonine at position 128 in its sequence this is a match. If, however, the annotated residue and the canonical residue didn’t match but another phosphorylate-able residue was found there instead, this was called a mismatch; for example, the Reactome annotation on Protein X indicates a phosphorylation on a threonine at location 362 whereas UniProt has a tyrosine at position 362. Finally, if the annotated residue was neither a match or a mismatch, the sequence was searched to find the closest exact match and the offset recorded. Once again using the example above, in Reactome another annotation on Protein X indicates a phosphorylation on a serine at position 495, whereas UniProt lists the closest serine to that location is actually at position 490; this would be counted as an offset of 5. It should be noted that if there was no match to the canonical residue all isoforms were then searched for the same criteria (match, mismatch and offset). Annotated modifications where the position is larger than the canonical sequence length, or where there is no position annotated at all, are considered to have not matched and no mismatch or offset is searched for.

Experimental data

We used phosphoproteomics data from 191 publications across 484 conditions obtained downloaded from the qPhos repository to assess the coverage of phosphoproteomics data by Reactome, PhosphoSitePlus, HPRD, BioGRID and SIGNOR. Data was downloaded from the qPhos download portal in January 2019. Each modification in qPhos was associated with a PubMed ID (PMID) and each PMID was tagged with metadata, such as the ‘sample’ (cell line) or condition the experiments were performed in. We subset the phosphosites by cell line and tissue type. Only phosphosites (and inherently PMID’s) tagged with a single ‘sample’ type (such as MCF-7 or HeLa) were analysed (in cases where there were multiple conditions it was impossible to discern which modification came from which sample).

Software and data availability

The source code and precompiled binaries are available from github (https://github.com/DavisLaboratory/Db_Compare). The software is compatible with Linux and MacOS operating systems and has been tested on both. The software requires Java version 8, Python 3 and R 3.6. Installation of these and dependencies is managed through the Anaconda package management system; a full list is available in the github repository along with installation instructions.

What proportion of the proteome is covered by different pathway databases?

The full extent of human phosphorylation is unknown and is an active field of research today. One study focused on profiling a single cell line identified approximately 75% of proteins to be phosphorylated (Sharma et al., 2014). Yet, many phosphoproteomics studies perform pathway level analysis at the protein level. This is due to the fact that many databases do not offer phosphosite level mapping. This motivated us to first determine the proteome coverage (and by proxy the coverage of proteins which may potentially be phosphorylated) by each database. This was done in order to provide context around their phosphoproteomic coverage as well as to assess their individual and combined potential for analyzing phosphoproteomic datasets. To do this, Reactome, WikiPathways, KEGG, HPRD, BioGRID, SIGNOR and PhosphoSitePlus were downloaded and their protein lists extracted. We normalised all protein UniProt identifiers to the January 2021 version of UniProt as explained in the methods. Because our focus is on resources that capture signalling reactions, throughout the analysis we used the subset of PhosphoSitePlus that contained kinase-substrate pairs (14,159 pairs) (henceforth referred to as PhosphoSitePlus (K-S) and PSP (K-S) in figures) and did not compare to phosphorylation events with no known kinase (231,645 events).

To characterize each database’s coverage of the proteome, a collection of reference knowledgebases and datasets were chosen for comparison (Fig. 1A). This set included qPhos, IMEx (Orchard et al., 2012), the Gene Ontology Cell Signalling Category (GO:0023052), the Gene Ontology Protein Kinase Activity Category (GO:0004672) (Ashburner et al., 2000) and the subset of UniProt annotated with the term ‘Phosphoprotein’ (KW-0597). The intersections between each reference dataset and each signalling database were then taken. Intersections with qPhos are expected to represent each knowledgebase’s coverage of phosphorylate-able proteins, while the intersections with IMEX are expected to represent each knowledgebase’s coverage of global interacting proteins, many of which have potential to be phosphorylated. We found BioGRID to have the highest coverage of each reference dataset and PhosphoSitePlus (K-S) had the lowest (barring the GO Protein Kinase Activity category where it had the 3rd highest coverage). Also of interest is the pattern of decreasing coverage of the all reference datasets across the signalling databases excepting the GO Protein Kinase Activity category. We see that BioGRID contains the most known human kinases, followed by SIGNOR, PhosphoSitePlus (K-S), WikiPathways, HPRD, Reactome and finally KEGG. Taken together these results point to a relative ranking for the analysis of data when mapping on the protein level, where you can expect BioGRID to have a high coverage of a given dataset, followed by Reactome, HPRD, KEGG, WikiPathways, SIGNOR and PhosphoSitePlus (K-S) but a differing coverage of the known human kinome. The same order was found when looking at each database’s coverage of SwissProt as well at the entirety of UniProt (Fig. S3, Supporting Information).

To explore the distinctions between each database’s proteome coverage we compared the database protein sets with each other using an UpSet plot (Conway, Lex & Gehlenborg, 2017). Overlaps of less than 50 were excluded, the entire UpSet plot can be found in Fig. S4, Supporting information. Examining the proteins common to all databases, we see the global intersection is just 5% of the total 19,816 unique proteins held across them (Fig. 1C), indicating a marginal global overlap. Individually, BioGRID consistently has the highest overlap with all other databases (Fig. 1B). Although the global intersection is small, the proportion of each database that is unique to itself ranges from only 2.5–20% (Fig. S4). We see that once again BioGRID consistently has the highest coverage of the other databases followed by Reactome in every database. WikiPathways and KEGG appear to have similar coverage of each other database. Finally, SIGNOR consistently has the 6th highest coverage (in all databases excepting PhosphoSitePlus (K-S)) and PhosphoSitePlus (K-S) the lowest coverage. The balance of proteins found in other databases vs. found only in a single database is an indication of potential information gain and novel insight when choosing which database and type to use. The set differences of each pair were also investigated to understand the scope of uniqueness of each reference dataset and signalling knowledgebase (Fig. S5, Supporting Information). Lastly, of note, PhosphoSitePlus (K-S) always maintains the lowest coverage of the other databases yet is the highest covered itself. This comparison is important as PhosphoSitePlus contains the most comprehensive known K-S network and each database’s relative coverage provides further evidence to its ability to connect phosphorylation signalling to a greater context.

As a potential consequence of the differing curation practices each database has adhered to, each database has tended to focus on a different area of biology. Of the number of proteins unique to each database, BioGRID had the most at 5,158, followed by Reactome with 393, KEGG and WikiPathways with 86 and 67 respectively, and then SIGNOR, HPRD and PhosphoSitePlus (K-S) all with less than 25 (Fig. 1C, coloured bars, Fig. S4 coloured bars). This alludes to the potential impact of database choice on conclusions drawn when performing pathway enrichment analysis. To explore this further, we performed an over enrichment analysis using the GO Molecular Function category on the proteins found to be unique to each database (Fig. S6, Supporting Information). We found that the proteins unique to Reactome were mostly involved in transporter activity. This could be because it is the only database in our comparison with abundant subcellular location annotations. BioGRID’s unique proteins were enriched for demethylase activity. WikiPathways was enriched for receptor activity. HPRD, SIGNOR, PhosphoSitePlus (K-S) and KEGG’s unique proteins did not produce any enrichment.

What proportion of the phospho-proteome is covered by different pathway databases?

We sought to address the question of the level of phosphoproteomic coverage of each database. None of the included databases (excepting PhosphoSitePlus (K-S)) have the ability to input a list of phosphosites and have them mapped on the phosphoproteomic level. To do this, all phosphorylation annotations were taken from each database where possible; neither KEGG nor WikiPathways were included in this section of the analysis as modification annotations in these databases are inaccessible.

In the first level of database characterization of coverage of the phosphoproteome, the databases of interest were compared to qPhos’ phosphorylation annotations as well as UniProt’s. qPhos contains 39,652 less phosphorylation annotations than the entirety of PhosphoSitePlus. However, qPhos was chosen as the authority on experimentally derived phosphorylations in this analysis, to avoid the self-comparison that would result from using a component of PhosphoSitePlus in evaluating its own utility. Pairwise Intersections were taken between each signalling knowledgebase and each reference dataset (Fig. 2A). Once again, we see a relatively consistent ranking of the databases, but in a differing order than the protein level. Here, HPRD has the highest coverage of qPhos’ phosphorylations, followed by PhosphoSitePlus (K-S), BioGRID, SIGNOR and then Reactome. A slightly different pattern occurs in the coverage of UniProt phosphorylations where the pattern follows the number of phosphosites in each database. Although HPRD has the highest number of phosphorylations annotated in the database category (with 43% of its phosphosites also recorded in UniProt) its coverage only represents 33% of all phosphosites recorded in UniProt. qPhos has almost 5 $\times$ more phosphorylation annotations than UniProt, and we see that HPRD has the highest overlap with 81% of its phosphorylation events found in qPhos, but still only representing 13% of qPhos’ total. Next is PhosphoSitePlus (K-S), SIGNOR and then BioGRID with 52%, 54% and 77% of their phosphosites respectively in common with qPhos. Yet, each database only covers a small proportion of qPhos at less than 3% for all three. Reactome has the smallest coverage of qPhos at 0.3% but this accounts for 46% of its own database. This means that approximately half of each of the latter 4 databases or more appear in global experimental datasets, but only a tiny proportion of the global experimental data can be interrogated through these knowledgebases. Taken together, this points to the disproportionate amount of data available and the need for tools to assign functional annotations. Added to the knowledge that only 20% of kinases are associated with 87% of currently annotated phosphosites (Needham et al., 2019), there is clear need to start to explore lesser known kinases and their contribution to cellular regulation.

Unlike on the proteomic level we see is no clear pattern in the overlap between databases (Fig. 2B). We see that the levels of uniqueness are much higher. The varying levels of overlap once more indicate the potential new information to be gained when mapping to a database. Keeping in mind that the subset of phosphosites used here from PhosphoSitePlus (K-S) all have a kinase or phosphatase linked to the feature. Though mapping events at the phosphorylation level is not yet supported by these databases, manual inspection of proteins identified in a study may reveal the presence of particular phosphorylation events of interest and their associated reactions, which may be then further explored. This analysis (Fig. 2) can be used as a guide to determining relative likelihood of finding a phosphorylation of interest in each database and the likely level of information that will be associated with it.

When looking at the inter-database coverage we see that the number of phosphosites found in all 5 databases is only 42, a mere 0.1% of the total 41,629 unique phosphorylation annotations across these databases (Fig. 2C) meaning each database covers a very different portion of the phosphoproteome. HPRD has 84% of its phosphosites unique to itself while PhosphoSitePlus (K-S) has 46%, followed closely by BioGRID at 63% Reactome at 37% and lastly SIGNOR at 13%. This is a clear indication of the vastness of the phosphoproteome and further the importance of database choice when it comes to analysing phosphoproteomic data. When looking at the number of unique phosphorylations belonging to each database, HPRD possesses the most at 26,423, PhosphoSitePlus (K-S) at 3,980, followed by BioGRID at 3,059, SIGNOR with 638 and Reactome with 599. The large number of unique phosphosites found in HPRD indicate its dedication to the curation of novel phosphosites (Goel et al., 2012). The set differences of each pair were also investigated on the phosphorylation level to understand the scope of uniqueness of each reference dataset and signalling knowledgebase (Fig. S7, Supporting Information).

Exploring the consistency between annotated phosphorylated residues and UniProt sequences

A major issue often found in phosphoproteomics is the incorrect assignment of substrates to kinases. This phenomenon is detailed by Humphrey et al. where they described that the correction of low-quality annotations is more difficult than an initial assignment due to the amount of evidence needed to refute prior knowledge. A major compounding factor to this problem is the lack of consistency in substrate phosphorylation annotations across databases. These mis-annotations can often be propagated by both human curators as well as text mining systems if left unchecked, contributing to the incorrect kinase assignments (Bachman, Gyori & Sorger, 2019). To explore this aspect of the problem we decided to test the consistency, or accuracy of the identity of annotated phosphorylated residues against canonical protein sequences found in UniProt. For example, if a database recorded a protein to have a phosphorylation of a serine on residue 456, we wanted to know if a serine would actually be found in the canonical (or isoform) protein sequence at position 456 (described in more detail in methods). Only phosphorylation sites annotated on the canonical protein sequence of a gene were included in Table 1. A match in an isoform was only searched for if a match was not found in the canonical sequence. Further, although a number of phosphorylations have been found on Histidine (Fuhs & Hunter, 2017), here we limited phosphorylate-able residues to Serine, Tyrosine and Threonine.

Mismatched sites can typically occur because their positions were annotated on post-translationally cleaved proteins, non-canonical isoforms, or on non-human orthologues (Bachman, Gyori & Sorger, 2019). It is likely that an offset of −1 is due to initial methionine cleavage. Offsets of larger size and differing direction are likely to be caused by other issues. The closest matching residue may not be the modified residue. In general, resolving inconsistencies like these will require individual investigation, however automated approaches have been developed that may address the most common sources (Bachman, Gyori & Sorger, 2019). Alternatively, allowing for small offsets (<2) may improve mapping by accounting for slight inconsistencies in annotations (Babur et al., 2018). We report all offsets here for completeness.

We found that the entirety of the PhosphoSitePlus database had the highest sequence consistency with UniProt canonical sequences at 99.7%. We found 260 cases where a match was found in an isoform rather than the canonical sequence, mostly owing to discrepancies with UniProt nomenclature (Table 1). In addition, we found 64 cases where the wrong phosphorylate-able residue was found, labelled a mismatch. Further, looking at PhosphoSitePlus’ K-S subnetwork, we see that it is almost the same at 99.6% consistent. We discovered only three phosphosites were found in isoforms as well as only 26 cases where the annotated residue was found next to the annotated location (i.e. offset by one residue). Next looking at qPhos we see a 99% percent consistency followed by SIGNOR with 98.5% consistency. So far each of these databases has shown exceptional curation practices in gathering phosphoproteomic data. Turning to HPRD we found overall phosphosite consistency to be 98.4% in relation to UniProt canonical sequences, followed by BioGRID with a phosphosite consistency of only 89.5%. Finally, we see Reactome which is 77.7% consistent. This lower consistency may be attributed to the manual curation practices of Reactome, or imprecise recording of isoforms in the literature curated for the database.

How well do databases cover experimental phosphoproteomic data?

Using phosphoproteomics to interrogate how cells sense and respond to their environment can lead to important insights into biological mechanisms behind these changes. Searching ‘phosphoproteomics’ in PubMed reveals over 3,500 results as well as depicting a steady rise in publications related to the search term over the last 20 years with 456 new publications in 2019 alone (Fig. S8, Supporting Information). As a result of the rise in interest and potential of phosphoproteomics, a number of databases have attempted to assemble the vast amounts of phosphoproteomic data being generated. As seen above, qPhos is one such database that provides curated phosphoproteomics datasets across 191 publications under a wide range of conditions. In this analysis we compared data from these cells with each knowledgebase to establish what proportion of the measurements from each cell type would map onto the knowledge contained in pathway databases. To pose the question another way, we sought to establish the coverage of pathway databases over phosphoproteomics experiments. Here we present results for the top six cell lines, ranked by the number of publications archived in qPhos (Table S3, Supporting Information). The full results are presented in Tables S4 and S5, Supporting Information.

Consistent with the coverage analysis we performed against the full proteome and the phosphoproteome, we see the same relative ordering in the databases in terms of coverage over experimental data. When compared at the protein level BioGRID has the highest coverage, followed by Reactome, HPRD, KEGG, WikiPathways, SIGNOR and lastly PhosphoSitePlus (K-S) (Fig. 3). Noting the order in which the cell lines are covered by each database, we see that coverage is highest for HELA cells across all databases, followed by HEK293, HeLa S3, Jurkat, MCF-7 and finally MCF-10. This order follows the number of proteins identified in each cell line in a decreasing order, which likely explains the pattern.

When comparing coverage of individual phosphorylation events, we see that HPRD has the highest coverage, followed by, BioGRID, PhosphoSitePlus (K-S), SIGNOR and finally Reactome (Fig. 4). In contrast to findings on the proteomic level, when looking at the phosphorylation level data we see a different pattern emerging in the cell line coverage by each database, where once again HeLa is the highest covered but Jurkat, HeLa S3, HEK293 and MCF-7 do not remain consistent and MCF-10 is the lowest covered (Fig. 4). Interestingly this pattern does not follow the number of phosphorylations per cell line. This may be due to the curation practices of each database, potentially highlighting underlying biases towards different common cell lines in experimental data from which pathway knowledge is drawn. These results also demonstrate that the proportion of measured phosphosites for which reaction level information is available is relatively small. For the most commonly measured cell line, HeLa, less than 1% of phosphosites map to annotations in Reactome and attach to known signalling pathways. This small number none-the-less corresponds to 22% of the phosphosites in Reactome. While HeLa phosphosite mapping to PhosphoSitePlus (K-S) is higher at around 4%, this knowledge doesn’t present the same level of mechanistic insight as is captured in Reactome, with functional information limited to K-S interactions. These results (Fig. 4) demonstrate that the vast majority of measured phosphosites lack pathway level annotation.

The contrast in coverage on these two levels by the same databases depict the consequences of generalizing phosphoproteomic level data to the protein level. Global coverage of qPhos can be compared in Figs. 1A and 2A and is reiterated in Figs. 3 and 4’s side-by-side comparison. Through this comparison we see that HPRD has a good coverage of phosphorylate-able proteins and a high coverage of phosphorylations with 56% of its proteins possessing a phosphorylation annotation. Conversely, BioGRID has a high coverage of phosphorylate-able proteins but many fewer phosphorylations when comparing a PPI of similar size (HPRD) with only 19% of its proteins recording a phosphorylation annotation. Reactome has a high level of phosphorylate-able proteins but a low number of phosphorylations translating to 5% of its proteins holding a phosphorylation annotation, while in PhosphoSitePlus (K-S) we see a low number of phosphorylate-able proteins but a high number of phosphorylations (96% of proteins with a phosphorylation annotation). SIGNOR, finally, appears to be on the lower end of proteins but possess a good number of phosphorylation annotations and has 33% of its proteins holding a phosphorylation annotation. This variation in the typical concentration of phosphorylations per protein emphasises the importance of database choice. It also highlights the different strengths of each database should for analyse phosphoproteomics data, as well as the biases that might result from database choice.