SNFM: A semi-supervised NMF algorithm for detecting biological functional modules

1.   Introduction

Proteins in cells seldom perform biological functions alone but form a larger molecular component with each other to exert specific functions ^[1,2]. With the development of high-throughput technologies such as mass spectrometry ^[3,4] and two-hybrid systems ^[5,6], numerous interactions among proteins have been acquired, resulting in large-scale protein-protein interaction (PPI) networks. Identifying the underlying protein functional modules from PPI networks is important to explore the biological process ^[7] in cells and to elucidate the disease-causing mechanisms ^[8]. Recent studies suggest that proteins densely interacting with each other are more likely to perform similar or the same biological functions ^[9,10], and various functional modules or protein complex detection algorithms have been proposed ^[11,12,13].

Typically, the PPI networks are modeled as graphs in which the proteins are the nodes and the interactions between the proteins are the edges ^[14]. A functional module can be regarded as sub-graphs with the nodes connected to each other more densely than others in the rest of the graph ^[15,16] with specific molecular functions. Thus, the goal of detecting functional modules from PPI networks can be determined using graph clustering method. For example, the cliques (fully linked sub-graphs) in a graph have been used for developing functional module detection algorithms, such as MC ^[17], LCMA ^[18], and CFinder ^[19]. Other proposed particular methods include core-attachment based methods, such as COACH ^[20]. However, these previous methods are only relied on the topological structure of the PPI networks. In addition, the current available protein interactions are rather limited and prone to be specious^[21]. Therefore, to improve the performance, Zhang et al. ^[22] proposed a generative model that takes the topological structure and Gene Ontology (GO) information together to detect functional modules. Georgii et al. ^[23] used gene expression data or phenotype data to build a weighted PPI network and then analyzed the functional modules in the PPI networks. Zhang et al. ^[24] developed an NMF-based functional module detection method named CoNMF that employs gene expression data and topological information simultaneously. To some extent, these methods reduce the drawback of those solely relied on the PPI network data, however, incorporating multiple data points would also induce more noises.

The databases on high-quality protein complexes, such as CORUM and PCDq ^[25], can be used for detect functional modules. Qi et al. ^[26] propose a semi-supervised functional module detection method named SCI-BN that extracts topological and biological features from known complexes and then uses them to train a Bayesian network for each sub-graph. Yu et al. ^[2] utilized the known protein complexes and the PPI network together to extract features and then used clique based algorithms to detect the functional modules. However, these two supervised approaches extract features from known complexes, which make them difficult to discover unknown functional modules.

To tackle these challenges, in this study, we propose a novel semi-supervised functional module detection algorithm(SNFM). Instead of extracting features from known complexes, we treated the protein pairs in one common complex as must-link constraints that can be viewed as prior information for our method. SNFM includes two main steps. The first step is to construct a similarity matrix of the proteins according to the adjacency matrix, in which we use a self-similarity manner to construct the similarity matrix as the featured matrix for the proteins. The second step is to factorize the similarity matrix instead of the adjacency matrix to obtain the protein module indicator matrix. Meanwhile, the protein pairs with must-link constraints are formulated as a graph regularization term. We will depict main components of the algorithm in detail in the Methods section.

2.   Methods

We first use the sparse subspace clustering method to construct the similarity matrix according to the protein adjacency matrix and next take the constructed similarity matrix as the protein feature matrix needed for the second step. The similarity matrix is constructed by taking into account the global information of the entire PPI network. Then, the similarity matrix is factorized into two non-negative matrices, and we propose an objective function based on NMF that takes protein pairs from complexes as must-link constraints. The similarity matrix obtained in the first step is decomposed to obtain the protein module indicator matrix, and the protein pairs of must-link constraints are formulated as a graph regularization term. We aim to improve the accuracy of the algorithm detecting the functional modules of protein interaction networks. The main steps of the proposed model are described in Algorithm 1.

2.1.   Subspace clustering

Given the PPI network G, G can be expressed as an adjacency matrix A $= (x_1, ..., x_N) \in \mathbb{R}^{ n\times n}$. The adjacency matrix is too sparse due to the limited protein-protein interactions. More importantly, the adjacency matrix only considers the relationship between the proteins and their neighbors and ignores the other topological information such as proteins with common neighbors that do not interact with each other. Each protein in the PPI network can be expressed as a linear combination of other proteins. Additionally, these n samples are from k subspaces $S_i$, i = 1, ..., k. The basic idea of sparse subspace clustering ^[27] is that a protein $x_i$, $x_i$ is represented as a linear combination of other proteins in the same subspace. Its corresponding objective function is as follows:

where Z$= [z_1, z_2, ..., z_n]$, Z$\in \mathbb{R}^{ n\times n}$ is the coefficient matrix. $Z_{ii} = 0$ is to avoid a protein expressed by itself. The larger the coefficient $Z_{ij}$ means the sample $x_i$ is more similarity to sample $x_j$; therefore, the obtained similarity matrix has superior subspace structure and robustness. However, in the specific application, there are noise data, or the sample is not clean, where E $\in \mathbb{R} ^{N \times N}$ is the error matrix.

2.2.   Supervised non-negative matrix factorization

Non-negative matrix factorization (NMF) ^[28] is a widely used matrix decomposition algorithm that can be used to detect overlapping protein modules in unsupervised manner. NMF decomposes the original matrix into two non-negative low-rank matrices, the product of which is the approximation of the original matrix.

However, most unsupervised module detection algorithms only consider the topological information of the PPI network. Moreover, due to the limited amount of PPI data obtained from high-throughput biological experiments ^[29] and the false positive links in PPI data, it is difficult to develop accurate functional modules by solely using PPI data. Currently, we can obtain information on high-quality human protein complexes from databases, such as CORUM ^[30]. Although in small scale, human protein complexes data is a high-quality data source of interaction network, which can be used as prior information. Therefore, we can design a semi-supervised functional module detection algorithm that combines the PPI data and the protein complexes data to compensate for the shortcomings of the PPI data.

The network similarity matrix A$\in \mathbb{R}_+^{n \times n}$ is used as the original input matrix. The main purpose of the NMF algorithm is to decompose the similarity matrix into the product of two non-negative low-rank matrices W $\in \mathbb{R}_+^{n \times k}$ and H$\in \mathbb{R}_+^{k \times n}$, where k$\ll n$. In this study, we use the Euclidean distance to measure the distance between Z and WH. Yang et al. ^[31] proposed a supervised community detection framework based on NMF and symmetric NMF. The objective function of the supervised model based on the NMF is as follows:

The matrix L is the Laplacian matrix, L = D - M. M is the must-link constraint matrix, and D is its corresponding diagonal matrix. $\beta$ is a positive parameter used to adjust the intensity of the prior information.

The column vector of the original matrix A is the weighted sum of all the column vectors in the left matrix W, and the weighted matrix is ​​the primary color of the column vector corresponding to the right matrix H; therefore, W is called the base matrix, and H is the coefficient matrix. The original matrix can be replaced by a coefficient matrix, and the features of the original matrix can be characterized and dimensionality reduced.

2.3.   Semi-supervised NMF Functional Module

In the original semi-supervised NMF model, the adjacency matrix only considers the information from the proteins and their neighboring proteins. Here, for SNFM, we use the similarity matrix constructed by sparse subspace clustering that takes into account the higher order neighbor information of the proteins, such as the protein neighbor's neighbor. We first use the subspace clustering method to construct the similarity matrix based on the adjacency matrix. Then, we performed the semi-supervised NMF algorithm with protein pairs from high-quality protein complexes as must links to improve the quality of the detected functional modules from the PPI network. Next we describe the algorithm of SNFM as follows:

where $Z = [z_1, z_2, ..., z_n]$, $Z \in \mathbb{R}^{ n\times n}$ is the coefficient matrix.

2.4.   Optimization

2.4.1.   Similarity matrix of PPI network

Subspace clustering is used to build the similarity of the PPI network that can be formulated as follows:

where A is the adjacency matrix of a PPI network and is a constant, and Z is the only variable. We can obtain the derivative of Z as follows:

Then using the KKT condition, we get the multiplicative iterative formula for Z as follows:

$\odot$ represents the product of the two matrix corresponding elements.

2.4.2.   Semi-supervised NMF

To minimize Eq. (2.3), we need to use the following matrix trace related knowledge: Tr (WH) = Tr (HW), $\|A\|_F^2 = Tr(AA^T)$, so Eq. (2.3) can be rewritten as follows:

Since W and H are both variables, Eq. (2.3) is not convex; therefore, to find the optimal W and H that can make Eq. (2.3) reach the minimum, we update one while keeping the other fixed as a constant. If we fix H as a constant, we can obtain a partial derivative of the objective function about W. The derivation process is reported as follows:

Similarly, if we fix W and treat H as a variable, the partial derivative of H can be calculated as follows:

Next, using the KKT condition, we obtain the updated rules of W and H as follows:

where D $\in \mathbb{R}_+^{n \times n}$ is a diagonal matrix of M, L = D - M is the Laplacian matrix of the matrix M, and $Tr(\cdot)$ is the trace of a matrix.

Matrices W and H are updated by the iterative iteration until the Eq. (2.3) converges or reaches the maximum number of iterations sets in advance. The time complexity is shown in Table 1.

3.   Results

To validate the performance of SNFM, we conducted experiments on three human-related PPI networks.

3.1.   Data and experimental setup

3.1.1.   Data sets

Three human-related PPI networks, namely the DIP, HPRD and STRING10, were used in our experiments. The DIP (Database of Interacting Proteins) ^[32] is a high-quality PPI database that integrates protein interactions from various data sources. It also contains interaction data calculated by algorithms from other reliable protein interaction databases with manually reviews. We extracted human-related PPI data from DIP, which includes a total of 4673 protein interactions and 2943 proteins. The HPRD (Human Protein Reference Database) ^[33] is a human-related database in which all of the information was manually extracted by biologists from the literatures. We extracted a total of 9453 proteins and 36888 protein interactions. The version of the STRING database ^[34] used in this study is STRING10, which is derived from several approaches, such as biological experiments, computational methods and literature curation. To obtain a high-quality PPI network, according to the introduction in ^[35], we only extracted the PPI data with a confidence score greater than 700 from the STRING10 database. Finally, we established a PPI network containing 14,380 proteins and 218,163 interactions.

In this experiment, we also used two human-related protein complex databases, CORUM and PCDq, as the validation data set. The CORUM (Comprehensive Resources for Mammalian Protein Complexes) ^[36] database provides data on protein complexes in human-reviewed mammalian organisms, all from the literature, excluding proteins obtained from high-throughput experimental data. We obtained the human related protein complexes from CORUM and removed the proteins that were not included in the human protein interaction data, so the number of protein complexes corresponding to different PPI data was different. PCDq is a human protein complex database. The protein complex data in this database is extracted from the integrated PPI network and then manually labeled for review and matched to complex data in the existing literature.

The corresponding CORUM protein complex in the DIP network contains 746 proteins, and PCDq contains 340 proteins. The corresponding CORUM protein complex in the HPRD network contains 1069 proteins, and PCDq contains 874 proteins. For STRING network, the corresponding CORUM protein complex contains 1134 proteins and PCDq contains 845 proteins.

3.1.2.   Evaluation metrics

We used Precision, Recall, and F1 to measure the quality of the detected protein functional modules. First, we defined the degree of coincidence between a protein complex $G_a$ and a detected protein functional module $D_b$ as follows:

where ${\left| G_a \right|}$ refers to the number of proteins in the protein complex $G_a$, ${\left| D_b \right|}$ refers to the number of proteins in the detected functional module $D_b$, and ${\left| G_a \cap D_b \right|}$ refers to the number of proteins appearing in both sets of $G_a$ and $D_b$. The larger the OL value, the more similar the two contrast sets are. In this study, if $OL(m, n) \ge \gamma$, we believe that the two protein sets m and n match successfully. According to the introduction in ^[37,38], we set $\gamma$ = 0.25 in this study. Next we defined the following two variables:

Subsequently, we defined the Precision, Recall, and F1 as follows:

where F1 is the harmonic mean of the precision and recall. The larger the indicator, the better the result of the detected functional module is.

For the human protein functional module, we first defined the neighborhood affinity score between the two protein sets as follows:

This value is used to measure the degree of overlap of the module p detected from the PPI network with a known protein complex g. If the neighborhood affinity values of p and g satisfy $NA(p, g) \ge \omega$, then we believe that p and g match. In this study, we set $\omega$ = 0.25. Suppose $m = |P|, n = |G|$, $t_{ij}$ represents the number of proteins that existed in the $i_{th}$ detected protein functional module and the $j_{th}$ protein complex simultaneously, and $N_j$ represents the number of proteins in the $j_{th}$ protein complex. Then the calculation formulas for Sn (Sensitivity), PPV (Positive Predictive), ACC (Accuracy) and MMR (Maximum Matching Ratio) are defined as follows:

The greater the ACC and MMR values, the greater the ability of the corresponding algorithm is to detect accurate functional modules from the PPI networks.

3.1.3.   Parameter analysis

The SNFM has only one parameter: $\beta$, which is used to control the smoothness between the topological information and the prior information. To investigate the influence of various parameters $\beta$ on performance, we set the value range of $\beta$ to 1, 2, 5, 10, 20, 50,100 and 1000. Taking the DIP database as an example, we used PCDq as the benchmark data set, and reported the changes of the corresponding ACC and MMR. As shown in Figure 1, we found that when the value of $\beta$ is approximately 1000, our proposed SNFM algorithm achieved better results. Thus, in all of the experiments in this study, we set $\beta$ = 1000.

3.1.4.   Baseline algorithms

We chose the two state-of-the-art algorithms for functional module detection: DCU ^[39] and OH-PIN ^[40], and four classical clustering algorithms, namely Markov clustering algorithm (MCL) ^[41], By finding the maximal clique, ^[42] has proposed an algorithm named CMC that can remove or merge clusters according to the interconnectivity. The NMF algorithm ^[28] is also used to detect functional modules. The k-means ^[43] algorithm performs clustering by initializing the cluster center as baseline methods.

3.2.   Experimental results

We obtained the related results for the six baseline methods and SNFM algorithm on three PPI networks. The corresponding Sn, PPV, ACC and MMR values of all comparison algorithms from three PPI networks (DIP, HPRD and STRING10) are shown in Table 2. For all comparative community detection algorithms, the algorithm proposed in this study SNFM includes a greater number of proteins in the functional modules and only slightly fewer than the DCU and k-means algorithm in the DIP network compared with other algorithms. The optimal ACC and MMR values were obtained by SNFM in all three protein interaction networks, and the number of detected modules that are matched protein complexes achieve the best results in all of the algorithms. The DCU, OH-PIN, MCL, and Myclique algorithms do not use priori information, and NMF and k-means use priori information. SNFM adopts the must-link method, and the effect is strong. According to ^[31], the reliability of the prior information is affirmed, indicating that the prior information is well utilized.

Additionally, we also analyzed the convergence speed of the SNFM objective function. In Figure 2, we show the trend of the value of the objective function as the number of iterations increases. We found that after a certain number of iterations (approximately 100 times), the value of the objective function tends to be stable; that is, we obtain the local optimal solution of the objective function.

Finally, using DIP, HPRD and STRING data sets, we obtained both the results using CORUM data as benchmark data(see Figure 3 (a), (b) and (c)) and those using the PCDq data as benchmark data (see Figures 3 (d), (e) and (f)). The results demonstrate that the SNFM algorithm proposed in this study has obtained the best results for precision(SNFM vs MCL for DIP database to improve 8%, recall(SNFM vs NMF for DIP database to improve 14%) and F1(SNFM vs NMF for DIP database to improve 12%).

3.3.   Functional module cases analyses

To investigate the biological insights of detected modules, we performed functional enrichment analysis using GO terms and listed the top 5 functional modules from DIP with significant biological processes (see Table 3). We found that there exist high ratio of proteins derived from same or similar protein complexes in the top 5 function modules. For example, in the M253, the proteins: PSMC2-6, PSMD3-4, PSMD6-7, PSMD10, PSMD12 and PSMD13, are from the multiprotein complex involved in the ATP-dependent degradation of ubiquitinated proteins. In addition, for the M263, which has 17 member proteins, we also found high ratio of proteins(e.g. MED6, 12-17, 19-20, 23-26) in a same protein complex. We further showed the corresponding subetaaph mapped from DIP network, which contains the 17 protein members of M263 and their first-order neighbors (see Figure 4). It showed that although the number of links between those protein members is rather high, there exist two obvious bipartite graph structures(Lee et al. ^[44]) in the subetaaph. We see that for the one part of the bipartite network, MED20, MED17, and CCDC belong to M263, but GATA1 and MED10 do not belong to M263. In the other part of bipartite graph, MED15 and MED6 belong to M263, but SREBF2 and SREBF1 do not belong to M263. This indicates that SNFM could detect the functional modules even with bipartite structures in a sparse network, such as DIP. This has been further validated by the additional enrichment analyses of pathway for the typical protein groups(e.g. MED17, MED20 and CCNC; MED6 and MED15, see Table 4) and GO terms for the whole M263(see Table 5).

SuperPaths is the unified GeneCards pathways that links to other pathways according to information extracted from the Kyoto Encyclopedia of Genes and Genomes (KEGG) and other databases. MED20, MED17 and CCDC are the first bipartite structures, MED15 and MED6 are the second bipartite structures, and CDK8 and MED24 are other members randomly selected in the module.

The more pathways shared by genes within the two map structures, the stronger the correlation is between them, which means the proteins should be clustered into the same module. Regulation of lipid metabolism by Peroxisome proliferator-activated receptor alpha (PPARalpha) and Metabolism are more special gene functions. This finding also illustrates the reliability of our proposed SNFM model for functional module detection.

We performed GO analysis on all proteins in M263. Table 5 lists the GO-terms including the following: biological process (BP), cellular component (CC), and molecular Function (MF). For example, the MF is GO:0001104(the RNA polymerase Ⅱ transcription cofactor activity), and the P-value is 2.16E-16, and the genes contained are MED19, MED15, MED26, MED17, MED24, MED14, MED13 and MED20, which account for 8/17 of all of the genes. Table 5 shows that MED17 and MED20 are highly correlated. Unfortunately, MED17 and MED20 have not interacted with each other in the DIP network, and most of the compared algorithms have not clustered them into the same module, except for k-means. However, the detected module by k-means, which contains MED17 and MED20, is too large (~1700 proteins). ^[45] has shown that larger modules have more diverse biological meanings.

4.   Conclusion

In this study, we propose a novel semi-supervised functional module detection model named SNFM, which uses the similarity matrix constructed by subspace clustering to not only consider the first-order neighbor information of proteins but also integrates global information from the PPI networks. The obtained similarity matrix is the feature matrix of the PPI network, and the use of limited prior information extracted from known reliable protein complexes allows for discovering the protein functional modules with significant biological meaning in the PPI networks. Through experiments on real human PPI networks (DIP, HPRD, and STRING10), we found that our proposed SNFM algorithm performs better in detecting protein functional modules. Because we currently have a high-quality protein complex with a small amount of data, our future work will consider designing a protein functional module detection algorithm that combines multiple biological datasets.

Acknowledgments

The work is partially supported by the National Key Research and Development Program (2017YFC1703506), the Fundamental Research Funds for the Central Universities(2017JBM020), the Special Programs of Traditional Chinese Medicine (201407001, JDZX2015170 and JDZX2015171), and the National Key Technology R&D Program (2013BAI02B01 and 2013BAI13B04).

Conflict of interest

We declare that we have no financial and personal relationships with other people or organizations that can inappropriately influence our work. We confirm that the manuscript has been read and approved by all named authors and that there are no other persons who satisfied the criteria for authorship but are not listed. We further confirm that the order of authors listed in the manuscript has been approved by all of us.

Supplementary

† These authors contributed equally to this work.