Predicting protein and pathway associations for understudied dark kinases using pattern-constrained knowledge graph embedding

Knowledge graph architecture

KGs are very similar to Heterogeneous Information Networks. An Information Network is a directed graph $G=\left(V,E\right)$, composed of vertices (also called nodes) and edges, with an associated node type mapping function ϕ:V → A and an edge type mapping function ψ:E → R (Shi et al., 2016). Each node v ∈ V belongs to one particular node type in the node type set $A:\varphi \left(v\right)\in A$, and each edge e ∈ E belongs to a particular edge type in the edge type set $R:\psi \left(e\right)\in R$. An Information Network is called a Heterogeneous Information Network if the sets A and R both contain more than one element, that is, there are multiple labels (types) for graph nodes and multiple labels (types) for edges. If sets A and R are singletons, the Information Network is called a Homogeneous Information Network (all nodes in the network are of the same type, and all edges are of the same type). While Heterogeneous Information Networks require that if two edges belong to the same edge type, the two edges share the same starting node type and ending node type, KGs do not. That is, the same edge (relation) type can be applied to different starting node types and different ending node types. This is also the case in the Resource Description Framework (RDF) (Cyganiak et al., 2014; W3C, 2014), a notation often used to represent KGs. RDF is based on the notion of triples of the form Subject-predicate-Object (Cyganiak et al., 2014) representing edges connecting nodes in the graph.

For example, in Fig. 1A, for a KG representing information about protein kinases, a node (entity) representing a dark protein kinase CDK 13 (cyclin dependent kinase 13) can be connected by an edge (relationship) labeled hasPathway to a node representing a pathway NeutrophilDegranulation, which represents the knowledge that CDK 13 participates in the pathway NeutrophilDegranulation. In such a KG, CDK 13 may be connected to other nodes using different labels (hasPathway), such as CDK13–hasMolelcularFunction—CyclinBinding, CDK13–hasBiologicalProcess—GranulocyteActivation and CDK13–hasCellularComponent—GolgiApparatus. Other protein kinases can additionally be included in the neighborhood context through shared nodes, such as the shared molecular function of CyclinBinding between dark kinase CDK13 and light kinase CDK6 (cyclin-dependent kinase 6). Here, edges have multiple labels, and destination nodes are of different types, which indicates that it is a heterogeneous KG.

Given a Heterogeneous Information Network $G=\left(V,E\right)$ (as defined above), the network’s schema is a directed graph, $S=\left(A,R\right)$, based on G’s node type mapping ϕ:V → A and its edge (relation) type mapping ψ:E → R.S is a directed graph defined over node types A, with edges as relations from R. Similarly, the schema of a KG represented in RDF is represented in RDFS (RDF Schema (W3C, 2014)). For example, given the CDK13 examples above, the schema would contain the edges Protein —participatesIn—Pathway, Protein–hasBiologicalProcess—BiologicalProcess, and Protein—isLocatedIn—CellularLocation. The schema, sometimes referred to as a meta-knowledge graph, or meta-graph, specifies constraints on using the edge labels to certain types of starting and ending nodes (subjects and objects in RDF). Also, given a KG schema, we can create a KG conforming to the schema containing many individuals (nodes).

Data sources and curation of datasets for knowledge graph generation

To construct our human kinase-centric KG, we gathered data from multiple publicly available curated resources (Fig. 1B). Nodes and edges within the graph were further integrated under more general node and relationship labels (also referred to as node types and edge types respectively) such as Disease (Figs. 1B and 2B). We further filtered the information included within our KG to avoid redundancy (described below). Fig. 1C shows basic statistics of different edges and their source and target nodes. The final KG contained 1,064,097 nodes of 11 types and 6, 279, 374 associations of 13 types.

The node types Protein and FunctionalDomain were populated with information on human kinase classification, functional domains, structure, and Protein–Protein Interaction (PPI). Information on these data types was extracted from the following databases: ProKinO (McSkimming et al., 2015), Dark Kinase Knowledgebase (darkkinome.org) (Berginski et al., 2021), UniProt (UniProt, 2021), PDBe-KB (Velankar et al., 2010), Pfam v33.1 (Mistry et al., 2021), and STRING v11 (Szklarczyk et al., 2019). Of note, included within our graph are numerous ortholog (non-human) Protein nodes collected from PPI databases. Protein Post-Translational Modifications (PTMs) associations were also included within our KG under node type PTM and were retrieved from iPTMnet v5 (Ross et al., 2017). The Chemical and Disease node types were populated from the Comparative Toxicogenomics Database (Davis et al., 2021) (protein-chemical, protein-disease, chemical-disease, and chemical-GO term associations) as well as cancer mutation data from the Catalogue of Somatic Mutations in Cancer (COSMIC) (Tate et al., 2019). The Genomics of Drug Sensitivity in Cancer (GDSC) (Yang et al., 2013) database was used for drug activity data. The Chemical node type additionally includes information on kinase-associated and ligand interaction motifs and ligand activity retrieved from Kinase-Ligand Interaction Fingerprints and Structures (KLIFS) (Kanev et al., 2021) and Pharos (Sheils et al., 2021), respectively. Gene ontology terms and pathway information regarding human kinases were retrieved from the Gene Ontology v2019_11 (Gene Ontology, 2021) (GOTerm node type) and Reactome v76 (Jassal et al., 2020) (Pathway node type), respectively.

Further filtering of data sources to reduce redundancy

The retrieved datasets were further curated prior to KG generation and population. All protein-pathway associations obtained from Reactome were split into manually curated associations (evidence = TAS) and predicted associations (evidence = IEA), only manually curated associations were included in our KG. Additionally, due to the hierarchical nature of the data taken from Reactome (i.e., Pathways are often organized into parent–child relationships), we eliminated high level pathways as they were deemed too general for our experiments. For example, we filtered out the pathways beneath “Disease” (R-HSA-1643685) but not beneath “Infectious disease” (R-HSA-5663205). This was done not only to reduce redundancy in the data but also to ensure the vector embeddings created from our KG for link prediction were not overfitted to highly connected nodes (often referred to as hubs). Similar filtering was done with Protein-GO term associations retrieved from Gene Ontology-GO terms with more than 5000 associations were removed. Additionally, Protein-GO terms were assigned separate node and edge types in the KG according to the MolecularFunction, CellularComponent, and BiologicalProcess terms.

Additional filtering to increase confidence in the data within our KG was performed using various strategies specific to the data source. Only PPIs (STRING) directly involving a kinase and with an experimental data score of more than 700 were included in our KG, and only PTMs (iPTMnet) with a confidence score of more than 1.0 were included in the graph. Additionally, Pfam entries with the protein kinase domains (“Pkinase” and “Pkinase_Tyr”) were removed.

Creation of curated kinase knowledge graph for kinase-pathway link prediction

The KG used in our link prediction experiments described here contains several node types such as Protein, GO Terms (BiologicalProcess, MolecularFunction, and CellularComponent), Disease, and other types. The nodes collected during a hypothetical random walk are shown in Fig. 2A. The hypothetical random walk is constrained by a schema organization that is shown in Fig. 2B the same schema is applied throughout the rest of the paper which indicates [Protein, any node not Protein and not Pathway, Protein, Pathway] (Jassal et al., 2020). In addition, the Protein type includes the understudied kinases (referred to as DarkKinases), well-studied kinases (referred to as LightKinases), and other proteins. This distinction is essential, as it gives us the ability to make predictions for a subset of proteins. The edge types (edges in the schema) are not labeled here, as we do not consider them in our method described here. Instead, we only rely on the types of source and target nodes. The loop edge returning to the Protein type indicates the Protein–Protein Interaction relation, which is included in our knowledge graph. Finally, Fig. 2C depicts the hypothetical subgraph in which the random walks are performed relating back to the nodes displayed in Fig. 2A.

Link prediction workflow

In the experiments and results presented in this paper, we used RegPattern2Vec. It is used as the first step in our link prediction process to produce vector representations for the nodes in the graph and formulate the link prediction as a classification problem. A ML model was trained using the combined vectors of existing pairs of nodes connected by an edge (link) of interest. The outline of the method is illustrated in Fig. 2D. We discuss each step of our link prediction method in the subsequent sections.

Regular pattern selection and its usage in random walks

Given a Heterogenous Information Network (as defined above), we have used regular patterns of the form H[^∧ T]+ H T for link prediction tasks, i.e., edges (links) between nodes of type H and T. In the pattern, H is a source node type and T is a target node type of the link of interest, respectively. Additionally, [^∧T] denotes any node type different from T (in regular expressions terminology it is known as the complement operator), and [^∧T]+ denotes a repetition (sequence) of one or more node types that are different from T. The overall pattern requires that any random walk must match the specified sequence of node types, in that a walk must begin with a node of type H then it must be followed by one or more nodes of type different than T, then a node of type H, which must then be immediately followed by a node of type T.

Simply put, the regular pattern is defined as a schema subgraph that should include only the relevant node types for a specific problem. An example regular pattern subgraph is shown in Fig. 2C. Unlike in the metapath2vec method (Dong, Chawla & Swami, 2017), the schema graph pattern represents multiple paths (walks) that can be used to predict missing links. Here, we aim to predict protein-pathway links and the random walks of interest are constrained by a regular pattern connecting node types and edge types most relevant to the prediction links. Here, when generating protein-pathway predictions, we have defined the regular pattern as Protein [^∧Pathway] + Protein Pathway, and each random walk instance must match it. Each walk starts from a Protein node. Selecting the next nodes on a walk is based on the existing graph nodes and matching the neighbor’s type in the regular pattern. When a random walk reaches a Pathway node, it follows the pattern in reverse order until a certain number of steps (nodes) in the walk are reached, or when it reaches a termination node, i.e., when there is no neighbor to match the pattern. We designed this specific aspect of the model due to the challenge of data sparsity for dark kinases. By allowing the reverse pattern to be followed we allow for additional relevant neighborhood nodes to be collected- when walking backwards once the pathway node is reached, the walk travels in a random direction until the hyperparameters are satisfied, capturing potentially relevant neighbors for our prediction task.

Biased random walk constrained by a regular pattern

As this work considered undirected graphs, enumerating a given graph’s path is computationally infeasible. The solution is to sample a subset of paths from all possible paths from a random distribution. The random walk constrained by a regular pattern is selected to generate walks of arbitrary length controlled by the “walk length” hyperparameter. Having the undirected heterogeneous network and a selected regular pattern, the random walk can be started from each instance of a starting node type in the pattern. As we want all the nodes to appear in our walks, iterating over them would be desirable. It is obvious that if we repeat the walk from each node (not all KG nodes may begin a walk, as constrained by the regular pattern), we will discover more walks as the node might link to multiple nodes of the same type. We call this hyperparameter the “number of walks” (referred to as NW in Results ‘Benchmarking RegPattern2Vec’s predictive performance with metapath2vec and other Graph Embedding Approaches’–Investigating the impact of hyperparameters on regpattern2vec’s predictive variability among replicate runs’). We will discuss how to choose the hyperparameter and the analysis of their impact in the Results ‘Benchmarking RegPattern2Vec’s predictive performance with metapath2vec and other Graph Embedding Approaches’–Investigating the impact of hyperparameters on regpattern2vec’s predictive variability among replicate runs’. The next step for each node is to select a node from the adjacent nodes based on the established regular pattern. This might result in multiple choices, and this is where randomization comes to play. RegPattern2Vec can utilize random distribution created by a user-defined function to generate the same probability for all the nodes or use an arbitrary distribution. To implement the random walk collection process, a regular pattern is converted to a Deterministic Finite Automaton (DFA), denoted by M. Each possible walk step is mapped to the DFA. The DFA M is used to check if transitions are allowed (an edge between two nodes); hence, a disallowed change gets a zero probability and is not used in the random walks.

On the other hand, in scale-free networks where the degree distribution follows the power law, some nodes may have a high degree of incoming/outgoing edges (known as hubs). Because such high degree nodes can dominate random walks and, consequently, representation learning, one popular way to address this issue is to bias the walks by the inverse of degrees of nodes (Grover & Leskovec, 2016), where the probability of choosing a node υⁱ⁺¹from υⁱ is calculated by normalizing the inverse of degrees of all neighbors of υⁱ. Although this approach lowers the probability of selecting high-degree nodes, it biases the random walk toward low-degree nodes, thereby capturing pertinent information related to low-degree “dark kinase” nodes. We used the formula below for node selection: ${\displaystyle P\left(v^{i+1}|v^i,M\right)=\left\{\begin{array}{c}\mathrm{g}\left(r^i\right)\frac{\frac{1}{\left|N_{v^{i+1}}\right|}}{\sum _{v\in N_{v^i}}\frac{1}{\left|N_v\right|}}\left(v^i,r^i,v^{i+1}\right)\text{is an edge in the graph and the}\hfill \\ \text{transition from v to v is allowed in M}\hfill \\ 0\left(v^i,r^i,v^{i+1}\right)\text{is an edge, but the transition is undefined}\hfill \\ 0\left(v^i,r^i,v^{i+1}\right)\text{does not exist in the graph}\hfill \end{array}\right.}$

In the formula, υⁱ denotes the current node and the candidate node for next step is υⁱ⁺¹ ∈ N_υⁱ.∣N_υ∣ denotes the degree of node υ,  and N_υⁱ is the set of all the neighbors of node ${\upsilon }^i.g\left(r^i\right)$ is the proportion of the edges of type r ⁱ among all edges of node υⁱ. Therefore, we randomly choose one edge (relation) type and then use the probability distribution by the inverse of node degrees to select the next node in the walk, but only among the nodes connected by the edge type chosen.

Vectorization of nodes and deep learning of semantic relationships between proteins and pathways

RegPattern2Vec uses a modified skip-gram model presented in (Dong, Chawla & Swami, 2017) to generate vector representations for the nodes of the KG. The random walks generate sequences of nodes, which resemble natural language sentences. The ML model simultaneously captures the local structure of the graph and the types of the nodes and encodes them as vector representations.

Link prediction as a classification problem

For each pair of protein-pathway associations, we combined their vector embeddings utilizing a widely used Hadamard product (Grover & Leskovec, 2016; Lin et al., 2015; Minervini et al., 2015). The resulting vector is used as features to train a Logistic Regression model. Training the model additionally requires generating negative examples for protein-pathway association pairs which are difficult as resources with true biological negatives, such as Negatome (Blohm et al., 2014), do not exist for our Pathway data type. To generate negative examples, we must first work under the closed-world assumption. This assumes that all information needed is provided, so for a statement to be true it must be explicitly stated to be so. Thus, a lack of a statement denotes that it is false (this contrasts with the open-world assumption in which statements can be true without them being explicitly known to be so). For our relation of interest (predicted link), we randomly select a head and tail node not connected by an edge in the graph. The embeddings of such nodes are then combined to produce negative examples. This process is repeated until the number of negative examples matches the positive examples.

Our logistic regression method utilizes the one-versus-rest strategy, in which a multi-class classification is split into one binary classification problem per class, meaning a classifier will be trained for each task, and therefore all the other data points in the split data are used as negative examples with the methods discussed above (under the closed-world assumption). Thus, the number of negative examples is dependent on the number of positive examples for each classification task. Although our dataset is large, this strategy is employed due to the inherent data imbalance between well-studied and dark kinases.

Capturing semantic relationships between proteins and pathways in deep learning

As previously mentioned in more detail, (see Materials & Methods ‘Vectorization of nodes and deep learning of semantic relationships between proteins and pathways’) RegPattern2Vec uses a modified skip-gram model presented in (Dong, Chawla & Swami, 2017) to generate vector representations for the nodes of the KG. Figure 3 shows the learned vector representation of nodes in the vector space using principal component analysis (PCA), a standard dimensionality reduction technique. For the protein-pathway predictions, we just consider the nodes to be of three types: “Protein”, “Pathway”, and “Others” when learning representation for the nodes. The separation of nodes in PCA shows our vector representation captures the node types and their network context. Several key differences exist between protein nodes and other nodes in our graph that may explain why they cluster separately from all other nodes. Protein nodes (including those annotated as dark/light kinases) are the default starting node in the schema, and many edge types describe proteins and connect them to other nodes in our graph (for example, a node describing a functional domain will be connected to the protein node it is describing and we additionally illustrate protein–protein interactions in our graph). This may allow protein nodes to occupy a unique space within the embedding reflected in the PCA.

Benchmarking RegPattern2Vec’s predictive performance with metapath2vec and other graph embedding approaches

To evaluate the accuracy of our method and note any improvements from metapath2vec, a test set was first generated by excluding 50% of the known protein-pathway associations from the training set. The training examples were generated using the process detailed in the Materials & Methods (Link prediction workflow–Link prediction as a classification problem’). To compare our method to metapath2vec, we produced AUC ROC curves for both models. Multiple curves were generated for the metapath2vec model as multiple meta-paths were used resulting in variability in the model’s performance (Fig. 4). Metapath 1 (‘Protein’, ‘FunctionalDomain’, ‘Protein’, ‘Pathway’, ‘Protein’) showed the best performance, based on 10-fold cross-validation, with an f1-score of 0.87 and AUC ROC of 0.94 for protein-pathway prediction. In contrast, RegPattern2Vec achieved an f1-score of 0.90 and AUC ROC of 0.96 on the same training and testing datasets (Fig. 4). We also compared RegPattern2Vec to other recently proposed embedding approaches for heterogenous graphs, namely BoxE (Abboud et al., 2020), MAGNN (Fu et al., 2020), and RHINE (Shi et al., 2020). These models were chosen mostly due to their ability to scale with larger knowledge graphs. Surprisingly, these models displayed poor performance in the task of predicting dark kinase pathway associations (Fig. S1), presumably because of the sparsity of knowledge related to dark kinases in the KG. RegPattern2Vec overcomes some of the challenges imposed by data sparsity by resampling a node as a starting point, allowing for constraint backward walking for a greater sampling of neighbor nodes, and use of inverse node degree as a consideration to select the next node in the walk.

Investigating the impact of hyperparameters on regpattern2vec’s predictive variability among replicate runs

Although RegPattern2Vec is better suited for the task of learning on sparse biological KGs, due to the nature of the random walks, there is an intrinsic level of variability with the predictions made by our model. Although the exact same sequence of nodes will not be sampled each time, there should be a larger context within the embedding space that identifies patterns useful for predictions for the dark kinases sampled. Therefore, overlaps in the Pathway-Protein predictions made by our model using different hyperparameters are expected. To demonstrate our model’s robustness, we measured the predictive consistency while modifying hyperparameters across three replicates (runs that all have the same hyperparameters) and compared these results. The percentage of protein-pathway predictions that overlap (within all replicates) was used as a metric for our model’s robustness. This analysis also revealed the variability in the overlap of the protein-pathway predictions on a per-kinase basis (for both light and dark kinases). This variability may be due to imbalances in the data available for each kinase. Notably, greater variability was seen amongst dark kinases, with the highest percentage overlap being 87.5% and the lowest percentage overlap being 0% (for the hyperparameter 40 NW with 95% confidence cutoff). A supplementary table (Table S2) provides further information indicating the overlap seen for each kinase. Additionally, we provided a supplementary figure (Fig. S2), displaying the confidence (0.6 to 1.0) in overlap for all hyperparameters tested.

One source of variability within our model is directly controlled by the hyperparameter “a number of walks (NW)”, which specifies the number of times a node should be sampled as the starting node during the random walk process (refer to Materials & Methods ‘Vectorization of nodes and deep learning of semantic relationships between proteins and pathways’). Repeating the walk from a specified node allows for a more complete characterization of neighbor nodes, often allowing for more node/edge types to be captured and considered in downstream analysis which is especially important when working with sparse data.

To find the optimal NW, the same starting node was used for generating paths in the random walk anywhere from 10 to 80 times. This process was repeated, resulting in a total of three replicates for each variation of the NW hyperparameter. Overlap between the protein-pathway predictions obtained through replicate hyperparameter conditions was then compared. From this analysis, we sought to identify the smallest value for the hyperparameter NW that can accurately capture local structure and node context within the graph to make consistent predictions (as determined by replicate overlap). This analysis revealed that changing the hyperparameter NW did not drastically change the number of overlapping protein-pathway predictions between replicate datasets (Fig. S3). Of note, the value of 40 for the hyperparameter NW (abbreviated 40 NW) generated the most consistent predictions among the NW parameters tested, with a mean of ∼22% overlapping pathways across three 0.7 cutoff replicates when averaging overlap over all kinases even when utilizing a range of confidence cutoffs from 0.7−0.95 (Fig. S3). As a result of this benchmark, we determined that 40 NW is the optimal hyperparameter for predicting protein-pathway relationships in the current KG, as it performs slightly better (as measured by replicate overlap) than the more computationally intensive hyperparameter values tested.

Hyperparameter optimization results in consistent replicate pathway predictions for 34 dark kinases across different replicate runs

After finding the optimal value for the hyperparameter NW (Fig. S3), we decided to investigate the overlap of all protein-pathway predictions produced by our previous analysis in which we test the hyperparameter NW (from 10–80 with a total of three replicate runs for each condition). The overlapping predictions present amongst highly variable walks suggest that these predictions are made based on nodes more often sampled within our network. After filtering for predictions that have over 0.95 confidence, only 95 dark kinase protein-pathway predictions emerged that were present in all replicates of all variations of NW tested (a total of 24 datasets) (Fig. 5A), resulting in representative protein-pathway predictions for 34 of the unique dark kinases in our dataset (Table S3). The top ten kinases (ranked according to the number of predictions that overlap) are shown in Fig. 5B. Of note, many of the protein-pathway predictions generated by our model were highly consistent (either in pathways with similar function or protein family involvement) on a per kinases basis (Table S3). For example, for the dark kinase VRK serine/threonine kinase 3 (VRK3), seven out of ten predictions were related to Toll Like Receptor (TLR) cascades, and five out of five protein-pathway predictions for the dark kinase protein kinase, membrane-associated tyrosine/threonine 1 (PKMYT1) were related to cell cycle control (Table S3). We further focused on a subset of dark kinases for path analysis described below.

Path analysis provides context for prediction associations between understudied PRKACB and DAG/IP3 signaling

The dark protein kinase cAMP-activated catalytic subunit beta (PRKACB) was amongst the kinases identified with the highest overlap of protein-pathway predictions (Fig. 5B). PRKACB, the beta-catalytic subunit of Protein Kinase A (PKA), appears to be involved in calcium regulation and modification, a known role for PKA (Barodia, Sophronea & Luthra, 2022; Kania et al., 2017; Ould Amer & Hebert-Chatelain, 2018; Palencia-Campos et al., 2020). Eight of the ten pathways identified were directly linked to calcium regulation (Table S3). To contextualize the predictions made by our model for PRKACB, we created a subgraph displaying common nodes traversed during the random walk process for the node representing the Reactome protein-pathway prediction R-HSA-1489509 (DAG/IP3Signaling). A portion of the subgraph highlighting the common nodes sampled between different paths are shown in Fig. 5C. The nodes found in all sampled paths are more likely to represent nodes that are important for providing the context that results in the shared protein-pathway prediction. The entirety of the subgraph (Data S1) displays all nodes with the addition of unconnected nodes that only appear in one path generated during the random walk process.

The highly sampled nodes of the subgraph shown in Fig. 5C, reveal that PRKACG (protein kinase cAMP-activated catalytic subunit gamma), PRKAR2A (protein kinase cAMP-dependent type II regulatory subunit alpha) (Manchev et al., 2014; Wei et al., 2021) and (calcium/calmodulin dependent protein kinase kinase 2) CAMKK2 (Najar et al., 2021) are nodes often sampled for the DAG/IP3Signaling prediction. The nodes being sampled multiple times demonstrate our model’s ability to leverage PPIs to infer pathway involvement. Additional nodes of interest that are often sampled for this prediction include PPIs between PRKACB (protein of interest) and both A-kinase anchor protein 28 (AKA28), and cGMP-specific 3′, 5′-cyclic phosphodiesterase (PDE5A). Both proteins have biologically relevant functions supporting the protein-pathway prediction for PRKACB’s involvement with DAG/IP3 signaling. AKA28 has been shown to bind to type II regulatory subunits of PKA (related to our dark kinase) and anchors/targets PKA to discrete locations within the cell (Kennedy & Scott, 2015; Kultgen et al., 2002; Omar & Scott, 2020; Sarma et al., 2010). PDE5A is a phosphodiesterase that regulates intracellular levels of cAMP and AMP (Peng et al., 2020). These observations support the predicted link/association between understudied PRKCB kinase and DAG/IP3 signaling pathway.

Predicted association of understudied CDK19 in TGF-beta receptor signaling: interpretability with path analysis and validation with PPI datasets

To further explore the biological relevance of our predictions for dark kinases, we compared our protein-pathway predictions to Reactome datasets generated using Protein Interaction and Proximity (PPI) data for dark kinases not included in our training data. PPI data was obtained from the Dark Kinase Knowledgebase (Berginski et al., 2021). Of note, many kinases in this dataset did not have sufficient protein-interactors to perform Reactome pathway enrichment analysis. For the 50 dark kinases with an adequate number of identified interactors, the overlap between the protein-pathway predictions generated by our link prediction and the Reactome enrichment generated by PPI data were compared. Figure 6A shows the top ten dark kinases ranked according to protein-pathway prediction overlap abundance and Table S4 provides the pathway information for all dark kinases with overlap seen. Among the dark kinases with the highest overlap from this comparison was cyclin-dependent kinase 19 (CDK19), a dark kinase belonging to the cyclin-dependent kinase superfamily (Malumbres, 2014).

After performing similar subgraph analysis to that previously described for PRKACB (Results,, ‘Hyperparameter optimization results in consistent replicate pathway predictions for 34 dark kinases across different replicate runs’), we were able to isolate nodes commonly traversed for the protein node of interest (CDK19) and the pathway node of interest TGF-betaReceptorSignaling (R-HSA-170834) (Fig. 6B). Again, only a subset of the paths is shown for ease of viewing. The entirety of the subgraph (Data S2) displays all nodes with the addition of unconnected nodes that only appear in one path generated during the random walk process. For the dark protein kinase CDK19 with the link prediction of TGF-beta receptor signaling, two “key” node hubs were identified including CCNC (cyclin-c) and CDK8 (cyclin dependent kinase 8). Recent total mRNA sequencing has demonstrated that loss of CCNC could activate the transforming growth factor (TGF)-beta signaling pathway (Tang et al., 2021). CDK8 has also been implicated in TGF-beta signaling through previous studies showing they drive Smad transcriptional action and turnover in TGF-beta pathways (Alarcón et al., 2009). Through interactions with several members of the mediator of RNA polymerase II transcription (MED) protein family, CDK19 can also indirectly exploit the relationship between CDK8 and TGF-beta receptor signaling. It appears that CDK19, and CDK8 have been shown to directly interact with not just each other but MED12, as well as CCNC, forming a complex collectively termed the “Mediator Complex” (Fant & Taatjes, 2019). This complex interacts with DNA-bound transcription factors and RNA polymerase II (Pol II) to activate and repress gene expression, with mediator subunit MED12 promoting TGF-beta signaling through both canonical regulation of transcription and non-genomic activity (Weber & Garabedian, 2018). Taken together this suggests that while direct connections contribute to the predictive power of the model the surrounding node context can also reinforce these connections and can contribute to predictions made.

Predicted association between understudied PSKH2 and cilium assembly

Protein-pathway predictions were additionally produced by our model for an understudied protein kinase PSKH2 (Byrne et al., 2022) which is a paralog of the Golgi-associated protein serine kinase H1 (PSKH1) (Brede et al., 2000). Currently, the Dark Kinase Knowledgebase (Berginski et al., 2021) lists only ten known protein interactors for PSKH2 and no pathway annotation for this kinase exists in the literature. We utilized this and another recently published dataset (Byrne et al., 2022) not included in training, for the Reactome enrichment analysis and independent validation of RegPattern2Vec predictions on PSKH2.

Comparing the overlap in protein-pathway predictions (40 NW dataset) with the Reactome enrichment, 18 pathways were shown to overlap (Table S5). The previously characterized subgraph visualization (Results, ‘Hyperparameter optimization results in consistent replicate pathway predictions for 34 dark kinases across different replicate runs’) was then performed, with five random paths chosen for visualization of the nodes commonly traversed for the protein node of interest (PSKH2) and the prediction node of interest (CiliumAssembly R-HSA- 5617833) (Fig. 7 and Data S3). Through this analysis one major “key” node hub was identified, UNC119B. unc-119 lipid binding chaperone B (UNC119B) is a protein that plays a key role in the localization of ciliary proteins by binding to N-myristoylated proteins and masking the hydrophobic lipid from hydrophilic cytosol, therefore facilitating trafficking to the primary cilium or the immunological synapse (Stroukov et al., 2019; Yelland et al., 2021). Many of the additional commonly traversed nodes belong to larger families of proteins known to play important roles in cilium assembly such as tubulins and dyneins. Thus, further characterization of these proteins can shed light on the role of PSKH2 in cilium assembly.

Overall, through our analysis of the paths utilized by the model as context for the link prediction process, we demonstrated that our modified random walk approach, RegPattern2Vec, can successfully sample multiple types of nodes and edges. These findings support our notion that our model preserves both local and latent structure through appropriate neighbor sampling and therefore provides the node context needed to make predictions. The predictions made by our model regarding dark kinases allow us to explore dark kinases of interest in a pathway context, giving us an idea of how these kinases may function and relate to other proteins.