Graph analysis and clustering of proteins linked with COVID-19

1.Introduction

Figure 1.

The Virion structure.

The current pandemic of COVID-19, a respiratory disease emerged in late 2019 has led to 1,051,635 confirmed cases and 56,985 fatalities in 208 countries with cases as of 4th April, 2020 [1]. The COVID-19 is induced by a new virus which causes severe acute respiratory syndrome – 2 (SARS-2) of the Coronaviridae family [2]. The history of coronavirus unfolds the episode of SARS-CoV in 2002 with 8000 confirmed cases and 10% fatality. Similarly, another event in Middle East respiratory syndrome coronavirus (MERS-CoV) in 2012 had 2500 confirmed cases and 36% fatality rate [3]. COVID-19 spreads with a higher fatality, and makes it hard to contain the disease and also escalates its pandemic potentiality [3]. From Latin, “corona” means crown. This remarkable pathogenic virus attack returns common flu to acute respiratory infections that can start long-term reduction in lung function, and proceed to death [4, 5, 6, 7]. Common flu spreads more quickly, and the fatality rate is low [8]. It is crucial to develop an understanding of how coronavirus proteins attack human proteins during infection, to devise therapeutic strategies to counteract COVID-19. This knowledge will be useful and can be applied to develop drugs and repurpose currently used ones. So far, no antiviral medications have been developed for COVID-19. The COVID-19 genome encodes with 16 non-structural proteins NSP1 to NSP16, four structural proteins E, M, N, and S, and eleven accessory factors, ORF3A, ORF3B, ORF6, ORF7A, ORF8, ORF9B, ORF9C, ORF10, ORF1AB, ORF1A, and ORF7B [9]. Table 1 shows the known COVID-19 protein names and description. Researchers are examining the presumed coronavirus proteins to develop drugs and vaccines in this pandemic. The Virion structure shows spikes projecting from the envelope that looks like a crown (Fig. 1) [10].

Protein-Protein Interactions (PPIs) present many challenges for the identification of drug-like molecules [11]. Both traditional as well as some innovative strategies like graph-based methods provided valuable tools for the discovery of PPI modulators and its potential of targeting PPIs for therapeutic intervention. We cannot understand the biological processes happening inside our body without extensive analysis of COVID-19 disease-human protein interactions [12]. The main objective of this work is to find protein clusters linked with COVID-19 disease-human interactions. Clustering proteins with the same biological functions help biomedical researchers to explore knowledge about the natural process happens in our body due to this disease. We propose the analysis of the structural and non-structural proteins as the seeds to construct a protein-protein interaction (PPI) network associated with COVID-19 [9]. In a PPI network, vertices or nodes represent proteins, and edges represent interactions. Graph cluster analysis with a blend of topological properties of the PPI network provides an appropriate biological knowledge for a promising tool to understand the biological function of the protein groups.

Further sections of this paper are in the following order: in Section 2, the authors did a literature study on the creation and investigation of PPI networks, Section 3 discusses how a PPI data can be modelled as a graph and the modelling algorithm which we have used to create the COVID-19 – human PPI data, Section 4 explains the graph clustering algorithm, Markov Cluster (MCL), Section 5 discloses the cluster validity measures, Section 6 describes the methodology of our analysis, Section 7 canvas the results and discussion, Section 8 covers possibilities, limitations, and future study, and finally, we conclude our study in Section 9.

2.Literature study

The literature study helps to accomplish a theoretical base for the research problem by seeking the past events of the particular subject. In this literature section, we cut across the knowledge about graph clustering and topological analysis of PPI networks. Topological means of PPI network analysis provide a platform for exploring complex diseases [13, 14, 15]. Network-based computational models to analyse COVID-19 disease – human protein interaction network – has already been applied by researchers in drug and therapeutics [16]. Studies used significant gene bio-signatures as the seeds to build the PPI network and analyse disease dynamics through topological analysis of the PPI network [17]. Several models for protein complex detection from PPI network has been applied in graph-based clustering techniques [18, 19].

Graph clustering using community detection algorithms follows a stochastic search method based on the seed vertex or edges [18, 20]. A heuristic graph clustering algorithm (HGCA) based on different topological properties has been proposed for protein complex detection [18]. The algorithm used a weighted degree for edges and vertices. A cluster description model is then constructed based on the candidate vertex and the cluster. Based on this description model, the HGCA starts with the seed vertex and produces communities in a greedy manner. A study with six community detection algorithms was used to analyse two biological networks and evaluated the resultant communities [20]. Among the six algorithms, Louvain cluster algorithm is found to be the fastest with an agglomerative approach to maximise the modularity. Density-based local neighbourhood search method and flow simulation method are widely used in protein complex detection and analysing PPI networks [19, 21, 22, 23, 24]. A demonstration of how MCL and Molecular Complex Detection algorithm (MCODE) has been done to identify patterns from PPI data related to Alzheimer’s disease [19]. A combination of MCL algorithm with Fruit Fly Optimization Algorithm (FOA), FOA-MCL has been used to find clusters formed by PPI network data of human immunodeficiency virus (HIV) [21]. An operon prediction model based on Markov clustering algorithm used some generic attribute information of genomes for graph [22]. The results show that the operon model has a better capability of operon prediction than classical operon prediction methods. The Random Walk Essential Proteins (RWEP) is a method that endorses random walks with restart that incorporate the topological and biological properties [23]. It has been applied to rule out protein essentiality in PPI networks.

Stochastic search clustering method is extensively used in PPI networks where the objects are considered as vertices. This method is time-consuming and, therefore, suitable for small networks. The flow-based method works on the principle of random walks and has a tendency to stay within clusters rather than between clusters. This method will not produce overlapping clusters. However, the process is time-consuming. Local neighbourhood density search method performs graph clustering by recognising seed proteins as individual clusters and then proceed greedily to add vertices.

Table 1 shows a summary of the reviewed graph clustering methods. Graph clustering and analysis using the clustering co-efficient as the validity measure provides a valuable tool for the partitioning of the PPI network [19, 25, 26]. However, the MCL algorithm is mainly designed for graphs and can be applied to biological applications [19, 27]. Novel algorithms for analysing PPI networks by combining MCL and optimisation techniques have also been developed [21].

3.Protein-protein interaction graph

A PPI network is a collection of protein interactions, often deposited in online databases [28, 29, 30]. Since proteins interact with each other and also carries information signals from one protein to another, understanding critical biological processes in the human body will be difficult without an extensive analysis of PPI [30]. Analysing PPI networks helps to mine data that assure one to create biological systems with new properties, and protein complexes for therapeutic purposes [16, 18, 31, 32, 33, 34].

Concerning topology, the PPI networks follow a small-world property and are scale-free networks [35, 36]. In small-world networks, it is possible to reach from a protein to any other protein in only a small number of steps. In scale-free networks, most proteins have a reduced number of interactions. Among the biological graph drawing algorithms, the force-directed layout model is very much flexible [37]. The algorithm can be used to calculate undirected graphs by using the information contained within the structure of the graph. There are several force-driven algorithms. Since the fundamental and aesthetic goal of optimisation of PPI graph drawing algorithms is the minimisation of crossings between edges, reducing the distance between edges and incorporating the properties such as small-world effect, scale-free network, the force-directed model is suitable for drawing a PPI graph [38, 39].

Generally, a PPI network can be created as an undirected, unweighted graph G=(V,E) where V is a group of proteins and E is a group of interactions between the proteins. The main data structure [40] used to store network representation is the adjacency matrix.

Let G be a graph such that V⁢(G)={v⁢1,v⁢2,…,v⁢n}, the adjacency matrix representation of G is a n×n matrix. If ar,c is the value in the matrix A, at row r and column c, then ar,c=1; if vr is adjacent to vc; otherwise, ar,c=0. Adjacency matrices require space of Θ⁢(|n|2).

The Kamada-Kawai [41] suggested that the number of edge crossings for a PPI graph is not a benchmark for a layout algorithm. This layout algorithm measures the total balance of the PPI graph as the square summation of the differences between the ideal distance and the actual distance for the entire vertices. The calculation is delineated in Eq. (1).

For some pair of vertices, i and j, where di⁢j is the standard distance amidst vertices, x is the set of coordinates and wi⁢j=di⁢j-α. The Kamada-Kawai approach preserves the total balance of the PPI graph, and deliver layouts with minimum edge crossings by approximation and minimisation of stress.

4.Graph clustering

Graph clustering is an augmenting area with a perspective to discover contemporary facts from complex data that can be pictured as a graph. The field of clustering has grown, and the number of clustering algorithms reported in biological applications is also high [19, 21, 42, 43, 44]. Graph clustering has two perspectives: intra-graph clustering and inter-graph clustering. Intra-graph clustering is the process of grouping objects within a single graph and inter-graph clustering method clusters between graphs. The intra-graph clustering method focus on both vector-based and graph-based. Vector-based clustering is of purely distance-based and graph-based clustering is established on the topological characteristics of the graph.

Graph-based clustering techniques for biological application explicitly uses theoretical graph terms to cluster the data that is represented as a graph. Local neighbourhood density search, flow-based simulation, and stochastic search are some of the graph-based clustering methods.

Figure 2 shows the perspective of graph clustering. Flow simulation method based on a random walk [45], or biological knowledge can be applied to PPI networks to uncover protein clusters. The MCL [27] is a flow simulation method based on a random walk. Moreover, the MCL is one of the most successful approach to cluster proteins in PPI networks [22, 46].

Figure 2.

Perspective of graph clustering.

5.Measures for cluster validation

Clustering is an unsupervised learning technique and gives distinct clustering results on the same data with distinct parameters. Clustering coefficient is one of the metrics available for some indication of the quality of the clusters [19, 25, 26, 47, 48, 49]. The clustering coefficient lies between 0.0 and 1.0. If the clustering coefficient tends to 1, the neighbourhood is fully connected, and graph possesses a maximal structure [50].

The dimension that represents the trend of a graph to be segregated into clusters is known as the clustering coefficient. A cluster is a subspace of vertices connected by edges.

Let v, be a vertex and let ev be the number of edges joining the jn neighbours of n. The clustering coefficient Cn of the vertex n, is delineated in the Eq. (5).

The clustering coefficient of CDi of a cluster Di is the average clustering coefficients of the entire proteins contained in Di. The clustering coefficient CD of the clusters D={D1,…,Dk} is defined in the Eq. (6).

For a biological network, clusters need to be validated using the domain knowledge to ensure the natural function of the objects. The most efficient way of biological validation is to test the gene ontology for the enrichment of classified clusters, and to check if proteins are functionally homogeneous.

Gene Ontology [51] is a database that consists of three categories of associations, namely, molecular function, cellular component, and biological process. The functions of proteins are labelled with a GO-term. A p-value to calibrate the biological annotation is computed with the hyper-geometric principle that is delineated in Eq. (7).

Figure 3.

The working process of the proposed model.

Figure 4.

COVID-19 disease – human PPI graph.

The proteins in the whole network and its linked GO-term are captured as the population (P) for the hypergeometric test. The proteins in the cluster with their matching GO-term is captured as the sample (S). This is to examine whether a specific GO-term enhances the cluster. The proteins in the population that are elucidated with a specific GO-term are the successes in the population (Ps⁢u). The proteins in the sample cluster that are elucidated with that specific GO-term are the successes in the sample (Ss⁢u). The threshold of p-value is set as 0.05.

6.Methods

The proposed model for finding proteins with related biological processes from a PPI network linked with COVID-19 disease consists of four phases. Figure 3 shows the working process for the proposed model.

6.1Data collection and preparation

The STRING is a database of predicted biomolecule interactions [52]. We used the database STRING to collect the PPI data. Also, the authors employed ‘Cytoscape’ automation to query the STRING database to retrieve network of proteins associated with COVID-19 disease [53].

Figure 5.

The adjacency matrix for created PPI graph.

The implicated COVID-19 disease proteins (“Nsp1”, “Nsp2”, “Nsp3”, “Nsp4”, “Nsp5”, “Nsp6”, “Nsp7”, “Nsp8”, “Nsp9”, “Nsp10”, “Nsp11”, “Nsp12”, “Nsp13”, “Nsp14”, “Nsp15”, “Nsp16”, “S”, “Orf3a”, “Orf3b”, “E”, “M”, “Orf6”, “Orf7a”, “Orf7b”, “Orf8”, “N”, “Orf9b”, “Orf9c”, “Orf10”) are given from ‘R’ tool using the package ‘RCy3’ [9, 54]. The RCy3 is an R package to link R language with Cytoscape. We restricted the query with the confidence score between 0.9 and 1.0 to get enough positive protein-protein interactions. The restricted disease protein query returned 57 proteins and 786 predicted interactions. The obtained COVID-19 disease network was then expanded with human host proteins. The authors kept the limit of host human protein interactors as 50. Here, we kept the confidence score as same as the disease protein-protein interaction network. The resultant final network contained 107 proteins and 1446 interactions.

Network visualisation of extensive protein-protein interaction data in a single frame is challenging. First, these networks tend to be large, typically consisting of hundreds of proteins with thousands of interactions between them. Cytoscape presents different layout algorithms. The authors created a force-directed spring embedded layout which uses the Kamada-Kawai algorithm to visualise the graph. The Cytoscape user interface allows us to drag nodes to override the automatic layout interactively. Figure 4 shows the constructed PPI graph.

6.3Clustering proteins

The MCL algorithm is established with random walks using Markov chain [19, 27]. The Markov chain is explained as, for a graph G, the data hoarded as a matrix M. Let r>0 be a number, the matrix derived after scaling each column of M with power coefficient r is called Γr⁢M, and Γr is called the inflation parameter with power coefficient r.

Write ∑r⁢qM for the summation of all the entries in column q of M raised to the power r. Formally, (Γr⁢M) is defined by the Eq. (8).

(8)

Γr⁢(Mp⁢q)=Mp⁢qr⁢/⁢∑r⁢q(M)

Every column q of a stochastic matrix M matches with the vertex q of the stochastic graph connected with M. The row entry p in column q matches with the probability of bustling from vertex q to vertex p. The clusters are formed when the matrix reaches a steady-state or all the values in a row become the same. The procedure for the algorithm is as follows:

Step1: Input: Adjacency matrix of the graph Step 2: add self-loops Step3: normalise the matrix to create Markov matrix Step 4: repeat steps 5 and 6 continuously to attain a convergence state Step 5: expansion by the e𝑡ℎ power of the Markov matrix Step 6: inflation of resulting matrix with parameter r>0 Step7: interpret the steady-state matrix to discover clusters

Figure 6.

The workflow of clustering proteins.

Figure 7.

Subgraphs of COVID-19 disease – human PPI graph.

Markov Cluster algorithm has been applied to the created COVID-19 disease – human PPI network. The workflow of clustering proteins is shown in the Fig. 6. During the clustering process, the expansion and inflation operations are performed repeatedly until a convergence state occurs. The expansion has been done by powering the matrix. The inflation has been done to prepare each element in the matrix. The overlapped cluster occurs when both clusters are graphically symmetric. It is not happening in this case.

The MCL returned eight clusters with an inflation parameter 1.8. Figure 7shows the obtained subgraphs for the clustered proteins.

7.Results and discussion

The COVID-19 disease – host human protein-protein interaction network with confidence level > 0.9 was downloaded from STRING database. The analysis has been done in RStudio environment on Intel Core-i5 5200U CPU, 64-bit 2.20 GHz processor and 8 GB of RAM. The R tool has over 6000 packages. The authors used igraph, mcl, and RCy3 packages for clustering and visualisation. The topological analysis of the PPI network has been done using igraph and tidyverse packages. The Table 2 shows the simple topological statistics of the created COVID-19 – human PPI network.

The structure and magnitude of the created PPI network are turning out to be 107 vertices and 1446 edges. The built PPI network should be further analysed to ensure that the network possesses scale-free property and small-world property. We analysed scale-free property of the constructed COVID-19 disease – human PPI network by power-law fit.

The degree distribution of biological networks approximates a power law: D⁢D⁢(d)∼p-d. This means that the probability or frequency of occurrence of a given degree in any vertex of the constructed PPI graph will be given as “p-d”, where p is the parameter characteristic of the PPI network and d is the numeric value of the degree. Figure 8 shows the frequency of degree distribution in the graph.

Figure 8.

Degree distribution of the PPI graph.

The proportion of variability in the degree distribution is computed on logarithmic value to fit the curve linearly and computed the R-square value. A negative exponential plot for the degree distribution of a network implies that the network follows power-law fit to predict the network’s scale-free property. Accordingly, Fig. 9 shows that the constructed PPI network’s degree distribution graphically fitted perfectly to a negative exponential plot. Besides, the R-squared value or coefficient of determination is reported as 0.906.

Figure 9.

Power law fit of the PPI graph.

In addition to graphical evaluation, we can adapt statistical evaluation also. The Kolmogorov-Smirnov test is carried mainly to evaluate scale-free networks that can be endorsed. According to this test, the power-law fit, p-value for the degree distribution of the PPI network is 0.9700409. The test result reveals that our constructed PPI network follows the scale-free property of biological networks. This strongly agrees that most vertices have a less degree and a few vertices have a higher degree.

Small-world properties of a network can be analysed with two topological properties: the average clustering coefficient or transitivity and the average path length, and its associated p-values. On the one hand, the transitivity of a vertex will be the fragment of probable edges between connected vertices that are literally represented in the network. On the other hand, the transitivity of all vertices can be averaged to get the clustering coefficient of the network. It will finally become the specific topological characteristic for the small world property of the network. High clustering coefficients exhibit the property of a small-world network, whereas small clustering coefficient indents the opposite. Here, the clustering coefficient or transitivity of the PPI network is 0.8683563. Another characteristic of a small-world network is the minimum average path length. Starting from a vertex, the minimum number of jumps or the minimum number of edges passes required to reach another vertex is known as the average path length. It is also known as the shortest path between the vertices. The average path length of our PPI network is 2.591233. Figure 10 shows the transitivity and average path length of our PPI network.

Figure 10.

Transitivity and average path length of the PPI network.

The clustering coefficient and average path length of the network alone does not imply that the constructed PPI network possesses small-world property. We need to calculate the p-values of the parameters, clustering coefficients and average path length. The p-value represents the probability of clustering coefficient and average path length score associated with a random network which is higher than our PPI network. The estimation of p-values can be done by calculating the clustering coefficient and average path length of the random networks and comparing them with our constructed PPI network. Therefore, the p-value can be estimated as the frequency with which the random networks overcome our constructed PPI network’s score. The accuracy of p-value depends on how large be the chosen random networks. In this example, the estimation of the p-value can be carried out with Barabasi game function in R, Sum (clustering coefficients > 0.8683563)/1000 = 0. The p-value associated with our PPI network’s average path can be drawn out from the same number of random networks. The sum (average path lengths > 2.591233)/1000 = 0.

Figure 11 shows that p-values associated with both parameters, transitivity and average path length have attained a score of 0 for thousands of random networks. The results of this evaluation strongly implicate that our constructed scale-free PPI network is also a small-world network.

Figure 11.

Graphical representation of a random small-world network.

A graph-based clustering technique using the MCL algorithm has been applied to the constructed COVID-19 disease – human PPI network. The MCL algorithm returned eight clusters without overlap. Among them, only seven clusters satisfied graph cluster validity through clustering coefficient between 0 and 1. Table 3 shows the results of MCL clusters.

Biological validation of the resultant clusters has been done using hypergeometric p-value test of GO-term. The gene ontology terms that are most common among the proteins that compose our eight clusters are found out. The test revealed that six clusters are involved in specific biological processes and satisfies the validity of p-value less than 0.05. The GO term finder returned the best feasible biological process associated with cluster (A): keratinization and annotated 38 proteins with p-value 3.50E-71. One protein has multiple mappings. The protein “LOR” has mappings to human protein identifiers “Q9UDR5” and “P23490”. Figure 12 shows the graphical representation of biologically investigated protein cluster (A).

Figure 12.

The graphical representation of protein cluster (A).

In cluster (B), 23 proteins are annotated to the regulation of cellular process with p-value 6.62E-05. There is one unmapped protein “PPP2R4”. Figure 13 shows the graphical representation of biologically investigated protein cluster (B).

Figure 13.

The graphical representation of protein cluster (B).

In cluster (C), the 22 proteins are annotated to the regulation of cell cycle with p-value 1.31E-27. Figure 14 shows the graphical representation of biologically investigated protein cluster (C).

Figure 14.

The graphical representation of protein cluster (C).

The whole proteins of cluster (D), cluster (E), and cluster (G) are annotated to mitotic cell cycle with p-value 1.66E-06, regulation of phosphoprotein phosphatase activity with p-value 1.15E-09, and G2/M transition of the mitotic cell cycle with p-value 3.03E-07, respectively.

The graphical representation of biologically investigated protein cluster (D), cluster (E), and cluster (G) is shown in Fig. 15. Meanwhile, cluster (F) and cluster (H) has no statistically significant results.

Figure 15.

The graphical representation of protein clusters (D), (E), and (G).

The authors created a COVID-19 disease – human PPI network and analysed the topological characteristics applying the MCL clustering to find protein clusters. Then the resultant clusters were analysed topologically and biologically. We relied on the evaluation measure and the clustering coefficient to validate the clusters topologically.

The clustering coefficient of the identified six protein clusters was either 1 or closest to 1. Hence, the sub-graphs possesses a maximal structure. The gene ontology term of proteins in each cluster is associated with a list of genes and results in a fine p-value < 0.05. The majority of proteins in the same cluster possess the same biological process. This implicates that each cluster is a protein complex with some biological significance. Consequently, the generated COVID-19 disease – human PPI network results in six statistically significant protein clusters with high clustering coefficient score. Figure 16 shows the statistics of proteins annotated to the gene ontology term from each cluster verses topologically significant clustered proteins. The number of proteins annotated to the gene ontology term for the biological process in each cluster, and the value of its clustering coefficient is pictured.

Figure 16.

Statistics of proteins annotated to the gene ontology term.

The topological analysis of the constructed COVID-19 disease – human PPI network provides evidence of the small world and scale-free properties of biological networks. Flow simulation-based graph clustering algorithm, MCL has been applied to the PPI network for further investigation and returned eight clusters. The clustering coefficient of MCL is equivalent to 0.8256, which implies that the clusters possess a maximal structure. The best feasible biological process associated with cluster (A) is keratinisation involving 38 proteins with p-value 3.50E-71. In total, 23 proteins of a cluster (B) are involved in the regulation of cellular process with p-value 6.62E-05, while the whole 22 proteins in a cluster (C) are involved in the regulation of cell cycle with p-value 1.31E-27. Meanwhile, the whole four proteins in the cluster (D), and cluster (E) are involved in mitotic cell cycle with p-value 1.66E-06, and regulation of phosphoprotein phosphatase activity with p-value 1.15E-09, respectively. The three proteins in the cluster (G) are reported to G2/M transition of the mitotic cell cycle with p-value 3.03E-07. However, cluster (F), and cluster (H) has no known significant biological terms. All the clusters were biologically validated using hypergeometric p-value test. Thus, six probable statistically significant protein clusters were identified.

8.Possibilities, limitations, and future study

The identified protein complexes spectacle knowledge about the biological functions in the disease – human interacting proteins. This gives an advantage to the COVID-19 researchers for therapeutic purposes. Also, the computational method using topological analysis would be a piece of valuable information for the researchers in other fields too. The authors have generated a medium-size PPI network. In this method, the clusters generated depend on two input parameters: expansion and inflation. These varying parameter values generate clusters with different size and structure. For an extensive PPI network, fine-tuning to generate clusters with good structure is a time-consuming process. This model can be enhanced by applying an automatic fine-tuning of the input parameter values according to the data set.

9.Conclusion

PPI network is a core aspect of the knowledge of biological organisms. Moreover, distinct physiological movements inside the human body are responsible for these interactions. The computational method using topological characters of the network helps in analysing the interactions between COVID-19 disease proteins and human proteins. In this paper, a computational analysis of the topological characteristics of the PPI network and graph clustering has been done. The clustering coefficient of the resultant clusters provides us with the information that the protein clusters possess a maximal structure. The presented model revealed that COVID-19 disease – human PPI contains groups of densely connected proteins involved in the same biological processes. Furthermore, these dense sub-graphs have a maximal structure with high clustering coefficients. The model finds valuable information about the biological process of the protein groups, which would be helpful for therapeutic researchers to understand the dynamics of COVID-19 disease. Consequently, it is expected that this study will provide a relevant contribution to the researchers in the field of biomedicine. In this model, the size and structure of resultant clusters depend on the varying tuning parameters, expansion, and inflation. For larger network, this is a time-consuming process. Our future work focuses on designing a more effective model by applying an optimisation algorithm to tune the parameters automatically in order to get the clusters from a large PPI network. Moreover, we will focus on incorporating additional biological information and different algorithms to discover protein complexes.