# Module Based Differential Coexpression Analysis Method for Type 2 Diabetes

Authors: Lin Yuan; Chun-Hou Zheng; Jun-Feng Xia; De-Shuang Huang
Journal: BioMed Research International
 (2015)

---

## Abstract

More and more studies have shown that many complex diseases are contributed jointly by alterations of numerous genes. Genes often coordinate together as a functional biological pathway or network and are highly correlated. Differential coexpression analysis, as a more comprehensive technique to the differential expression analysis, was raised to research gene regulatory networks and biological pathways of phenotypic changes through measuring gene correlation changes between disease and normal conditions. In this paper, we propose a gene differential coexpression analysis algorithm in the level of gene sets and apply the algorithm to a publicly available type 2 diabetes (T2D) expression dataset. Firstly, we calculate coexpression biweight midcorrelation coefficients between all gene pairs. Then, we select informative correlation pairs using the “differential coexpression threshold” strategy. Finally, we identify the differential coexpression gene modules using maximum clique concept andk-clique algorithm. We apply the proposed differential coexpression analysis method on simulated data and T2D data. Two differential coexpression gene modules about T2D were detected, which should be useful for exploring the biological function of the related genes.


---

## Body

## 1. Introduction

DNA microarray has been widely used as measurement tools in gene expression data analysis [1–4]. Gene expression profiling data from DNA microarray can detect the expression levels of thousands of genes simultaneously, providing an effective way for mining disease-related genes and revealing information of the regulatory networks and biological pathways of genes. Currently, the analysis of gene expression data can be divided into three levels: first, analysis of the expression level of individual genes, determining its function based on gene expression level changes under different experimental conditions: for example, the tumor type specific genes are identified according to the significance of difference in gene expression using the statistical hypothesis testing analysis method; second, study of gene interaction and coregulation through the combination of genes and grouping; and, third, an attempt to deduce the potential gene regulatory networks mechanism and explain the observed gene expression data.Among the microarray data analysis methods, gene differential expression analysis is one of the most widely used types of analysis for disease research. Gene differential expression analysis method selects differentially expressed genes according to expression change value of a single gene. In fact, gene expression value change between normal samples and disease samples can be used to present the possibility of the relation between gene and disease. However, the traditional pathogenicity genes selection methods based on gene expression data treat each gene individually and interaction between them is not considered. Actually, genes and their protein products do not perform their functions in isolation [5, 6], but in cooperation. Functional changes such as alteration in tumor cell growth process, energy metabolism, and immune activity were accompanied with coexpression changes. Differentially expressed genes selection methods often focus only on the size of the single genes and the relationship of individual genes and disease, ignoring a plurality of pathogenic genes of the complex disease as a gene module with disease related, as well as within the module gene [7].Differential coexpression analysis, as a more comprehensive technique to the differential expression analysis, was raised to research gene regulatory networks and biological pathways of phenotypic changes through measure gene correlation changes between disease and normal conditions. Differential coexpression genes are defined as genes whose correlated expression pattern differs between classes [8]. The gene coexpression changes between different conditions indicate gene regulatory pathways and networks associated with disease. In gene differential coexpression analysis, a pair of gene expression datasets under disease and normal conditions is transformed to a pair of coexpression matrix in which links represent transcriptionally correlated gene pairs [5]. Until now, methods for differential coexpression analysis of gene expression data have been extensively researched, and multiple algorithms have been developed and tested [9–12]. In those gene differential coexpression analysis methods, the most common choice of similarity measurement is Pearson’s correlation coefficients. However, Pearson’s correlation is sensitive to outliers. So biweight midcorrelation (bicor) is considered to be a good alternative to Pearson’s correlation since it is more robust to outliers [13].In biomedical research, many complex diseases are contributed jointly by alterations of numerous genes; they often coordinate together as a functional biological pathway or network and are highly correlated. With recent interest of gene differential coexpression analysis in the gene network or module, gene module analysis has emerged as a novel holistic approach for microarray analysis. Somewhat large units, made up of genes, are more densely connected to each other than to the rest of the network, are often referred to as modules, and have been considered to be the essential structural units of real gene networks. There exists overlap among gene modules in large real networks.Until now, there are many methods to find gene modules. For example, Butte and Kohane [14] proposed a systems-based approach called Entropy Minimization and Boolean Parsimony (EMBP) that identifies, directly from gene expression data, modules of genes that are jointly associated with disease. Kostka and Spang [15] used additive model to find differential coexpression gene modules. Prieto et al. [16] used altered expression based on improved additive model, optimal residual ratio, and minimumF-distribution to find differential coexpression gene modules. However, the microarray data contains a large number of genes; those methods need to search all gene expression data resulting in a large amount of computation; the process is very time-consuming even using optimized search algorithm.The maximum clique analysis can avoid exhaustive search and quickly find maximum gene module with biological significance. The maximum clique problem (MCP) is a classical combinatorial optimization problem in graph theory. In 1957, Ross and Harary [17] first proposed the deterministic algorithm to solve the maximum clique problem. Since then some researchers had presented a variety of algorithms to solve this problem. The maximum clique problem is widely used in different areas, such as signal transmission, computer vision, and biological research. In this study, a gene coexpression network can be treated as a graph; gene is represented by vertex and coexpression relationship is represented by edge. We will usek-clique algorithm [18], which is an effective and deterministic method for uniquely identifying overlapping modules in large real networks. We first show some basic definitions.k-cliques, the central objects ofk-clique algorithm investigation, are defined as complete (fully connected) subgraphs of k vertices.k-clique adjacency is as follows: twok-cliques are adjacent if they share some vertices.k-clique chain is as follows: a subgraph, which is the union of a sequence of adjacentk-cliques. We usek-clique algorithm to find gene cliques, and maximum clique concept is used to quickly find large gene modules which are made ofk-clique chain. For the sake of convenience, we use the terms graph and community or network interchangeably, the former stressing the mathematical concept and the latter the application.In this paper, we proposed a new approach for gene differential coexpression analysis in gene modules level based on combining biweight midcorrelation, differential coexpression threshold strategy, and maximum clique concept andk-clique analysis. Biweight midcorrelation measures the coexpression relationship between genes and thek-clique analysis with maximum clique concept quickly finds maximum disease-related module with biological significance. We use the approach to further investigate the gene module in order to gain insight into coexpression relationship between genes. The algorithm can find differential coexpression disease genes modules and global coexpression patterns are determined for type 2 diabetes expression dataset. As far as we know, no one has done this experiment.The rest of the paper is organized as follows. Section2 describes the methods proposed in this study. The biweight midcorrelation coefficients, “gene differential coexpression threshold” strategy, and threshold selection strategy are first presented, and the algorithm ofk-clique is consequently given. Section 3 presents the experiment on simulated data and type 2 diabetes (T2D) in rats dataset. Section 4 concludes the paper and outlines directions of future work.


## 2. Methods

### 2.1. Biweight Midcorrelation for Differential Coexpression

Differential coexpression analysis usually requires the definition of “distance” or “similarity” between measured datasets, the most common choice being Pearson’s correlation coefficients. However, Pearson’s correlation coefficient is sensitive to outliers [13]. Biweight midcorrelation is considered to be a good alternative to Pearson’s correlation since it is more robust to outliers. Example of a gene expression matrix is as follows:(1)Gene1Gene2⋯GenepZ1X11X12⋯X1pZ2X21X22⋯X2p⋮⋮⋮⋱⋮ZnXn1Xn2⋯Xnp.For each sampleZi, we measure expression levels for a set of genes, so Xij is the measurement of the expression level of the jth gene for the ith sample, where j=1,…,p. The xth column vector of matrix represents gene expression profile of gene X. In order to define the biweight midcorrelation (bicor) [13] of two numeric vectors x=(x1,…,xm) and y=(y1,…,ym), we first define ui, vi with i=1,…,m:(2)ui=xi-medx9madx,vi=yi-medy9mady,where med(x) is the median of vector x, mad(x) is the median absolute deviation of vector x, mad(x) is the median of new numeric vector in which each number is absolute difference between original vector value and med(x); this leads us to the definition of mad(x) and weight wi for xi, which are(3)madx=medxi-medX,wix=1-ui22I1-ui,where the indicator I(1-ui) takes 1 if 1-ui>0 and 0 otherwise. Thus, the weight wi(x) is close to 1 if xi is close to med(x), approaches 0 when xi differs from by nearly  9mad(x), and is 0 if xi differs from med(x) by more than 9mad(x). An analogous weight wi(x) can be defined for yi. Given the weights, we can define biweight midcorrelation of x and y as(4)bicorx,y=∑i=1mxi-medxwixyi-medywiy·∑j=1mxj-medxwjx2hhhih·∑k=1myk-medywky2-1.It should be noted that the equations of biweight midcorrelation do not involve an explicit identification of outliers, and all elements whose weightwi=0 can be considered outliers. The user can also set up the maximum allowed proportion of outliers using the argument “maxPOutliers”; the “maxPOutliers” is interpreted as the maximum proportion of low and high outliers separately. For the value of bicor from −1 to 1, −1 represents the maximum negative correlation and 1 represents the maximum positive correlation. Zero represents irrelevant correlation.


### 2.2. The “Differential Coexpression Threshold” Strategy

We used biweight midcorrelation to measure every pair of genes in the gene expression dataset and get a gene coexpression matrix. The gene coexpression matrix is a square and symmetric matrixP whose rows and columns correspond to the genes and whose element Pij denotes the coexpression relationship between genes i and j. In this paper, we use AGN which represents gene coexpression adjacency matrix in normal conditions and AGD which represents gene coexpression adjacency matrix in disease condition. To find differential coexpression gene modules which are coexpressed in normal condition and not related to disease condition, we set two thresholds T1 for adjacency matrix AGN in normal condition and T2 for adjacency matrix AGD in disease condition. AGN(i,j) is set to 1 if value of AGN(i,j) is greater than or equal to T1; otherwise, AGN(i,j) is set to 0 and AGD(i,j) is set to 1 if value of AGD(i,j) is less than or equal to T2; otherwise, AGD(i,j) is set to 0. We integrated AGN and AGD into a matrix AG after we had intersection of the corresponding elements of AGN and AGD. AG(i,j)=1 means coexpression value of gene i and gene j in AGN is greater than or equal to T1, and coexpression value of genes i and j in AGD is less than or equal to T2. AGN(i,j) also can be set to 1 if value of AGN(i,j) is less than or equal to T1 and AGD(i,j) is set to 1 if value of AGD(i,j) is greater than or equal to T2. The method is shown in (5).With the above mentioned strategy, we also setAG(i,j)=1 if the absolute value of AGN(i,j) subtracting AGD(i,j) is greater than or equal to T3 and the absolute value of AGN(i,j) is greater than or equal to the absolute value of AGD(i,j) simultaneously. This is a special type of coexpression change. In reality, coexpression reversal probably has biological significance. The coexpression reversal between normal condition and disease condition has advantage in disease. For example, the coexpression ofp53 andKlf4 recently reported that the positive or negative correlation between these two genes determines the outcome of DNA damage, DNA repair, or apoptosis [19]. We believe that our attention to this special coexpression change will help to explore subtle mechanisms involved in genes transcriptional regulation. We excavated maximum cliques which have biological significance from AG adjacency matrix to further investigate gene regulatory networks. Consider the following:(5)ifAGNi,j≥T1,thenAGNi,j=1,elseAGNi,j=0;ifAGDi,j≤T2,thenAGDi,j=1,elseAGDi,j=0;AGi,j=AGNi,j∩AGDi,jifAGNi,j≤T1,thenAGNi,j=1,elseAGNi,j=0;ifAGDi,j≥T2,thenAGDi,j=1,elseAGDi,j=0;AGi,j=AGNi,j∩AGDi,j.


### 2.3. The Threshold Selection Strategy

The two real value adjacency matrixes are transformed into a binary matrix which contains two elements 0 and 1 only. Choosing different thresholds will lead to different results; too largeT1 threshold or too small T2 threshold will lead to small link number, low density clique, and lost biological significance cliques. On the other hand, too small T1 or too large T2 will lead to many overlapping cliques. They are not helpful for finding biological significance differential coexpression gene disease-related modules. In fact, how to choose a reasonable threshold in conversion process is a problem which needs to be further studied. Generally, the selection of the threshold can be based on the proportion of outliers in the figure or the density of graph. The outlier is the point which is not connected to any edges. The density is defined as the ratio of number of edges to the maximum possible number of edges in the graph. The density of clique is 1.For gene expression data analysis, closely linked functional module is not the strict sense of maximum clique due to the lack of certain section. In this paper, we use density to measure approximation degree of functional module with gene differential coexpression clique, which may be having more biological significance.


### 2.4. The Maximum Clique Concept andk-Clique Algorithm

Graph theoretical concepts are useful for the description and analysis of interactions and relationships in biological systems. In gene coexpression graph, gene is represented by vertex and coexpression relationship by edge.G=(V,E) is an arbitrary undirected and weighted graph unless otherwise specified in graph theoretical concepts. V={1,2,…,n} is the vertex set of G, and E is the edge set of G. For each vertex i∈V, a positive weight wi is associated with i. AG=aijn×n is the adjacency matrix of G, where aij=1 if (i,j)∈E is an edge of G, and aij=0 if (i,j)∉E. Genes and relationship between genes are represented by vertex and edge, respectively.A graphG=(V,E) is complete if all its vertices are pairwise adjacent; that is, for all i,j∈V, (i,j)∈E. A clique C is a subset of V such that G(C) is complete. The maximum clique problem asks for a clique of maximum weight. An independent set (stable set and vertex packing) is a subset of V, whose elements are pairwise nonadjacent. The maximum independent set problem asks for an independent set of maximum cardinality. The size of a maximum independent set is the stability number of G (denoted by α(G)). The maximum weight independent set problem asks for an independent set of maximum weight. A maximum clique means a clique which is a subset of the nodes in V in which every pair of nodes in the subset is joined by an edge and is not a proper subset of any other cliques [20].In application, the identification of maximal cliques is often of limited interest since the requirement of complete connectivity is so restrictive. When dealing with imperfect systems or with experimental data, we may need to consider more general notions of cohesive subgroups. In this paper, we consider different notions of cohesive subgroups that includen-cliques, k-plexes, and λ-sets [18]. It is well known that the nodes of large real networks have a power law degree distribution [21]. Most real networks typically contain parts in which the nodes (units) are more highly connected to each other compared to the rest of the network. The sets of such nodes are usually called clusters, communities, cohesive groups, or modules [22–26], which have no widely accepted unique definition. The basic observation on which our modules definition relies is that a typical gene differential coexpression module consists of several complete (fully connected) subcliques that tend to share many of their nodes. To find meaningful communities, several basic requirements should be satisfied: it cannot be too restrictive, should be based on the density of links, is required to be local, should not yield any cut-node or cut link (whose removal would disjoin the community), and, of course, should allow overlaps. We employ the community definition specified above because none of the others in the literature satisfy all these requirements simultaneously [27–29].k-clique algorithm for detecting gene differential coexpression modules in a network has been published in the paper [26].k-clique algorithm is also named clique percolation method. The existing divisive and agglomerative methods recently used for large real networks have some disadvantages. Divisive methods cut the network into smaller and smaller pieces; each node is forced to remain in only one community and be separated from its other communities, most of which then necessarily fall apart and disappear [27, 30]. The agglomerative [31] method has the same problem. Thek-clique algorithm has demonstrated the advantages over the divisive method and agglomerative method. In the algorithm, although the numerical determination of the full set ofk-clique communities is a polynomial problem, the algorithm is exponential and significantly more efficient for the graphs corresponding to actual data. Thek-clique algorithm first locates all cliques (maximal complete subgraphs) of the network and then identifies the communities by carrying out a standard component analysis of the clique-clique overlap matrix [28]. Thek-clique algorithm uses the threshold probability d(k) (critical point) ofk-clique percolation to find all maximal complete subgraphs. The critical point is shown in (6), where N is the number of genes or vertex of graph:(6)dk=1k-1N1/k-1.Thek-clique algorithm gives two plausible choices to measure the size of the largestk-clique percolation cluster in (7) and (8). The most natural one, which we denote by N∗, is the number of vertices belonging to this cluster. ϕ is an order parameter associated with this choice as the relative size of that cluster:(7)ϕ=N∗N.The other choice is the numberL∗ ofk-cliques of the largestk-clique percolation cluster. The associated order parameter is again the relative size of this cluster:(8)φ=L∗L,where L denotes the total number ofk-cliques in the graph. L can be estimated as(9)L≈Nkdkk-1/2≈Nkk!dkk-1/2.In this paper, we use the biweight midcorrelation for constructing binary networks. Two-condition coexpression adjacency networks can always be transformed into a binary one by ignoring any directionality in the links and keeping only those stronger than a threshold weight. Then, the concept of maximum clique andk-clique algorithm were used to find gene differential coexpression modules. We named the proposed method “BMKC” (biweight midcorrelation andk-clique algorithm) method. Changing the threshold is like changing the resolution with which the community structure is investigated: by increasing, the communities start to shrink and fall apart. A very similar effect can be observed by changing the value of k as well: increasing k makes the communities smaller and more disintegrated but, at the same time, also more cohesive. More details aboutk-clique algorithm can be found in [28, 32].


## 2.1. Biweight Midcorrelation for Differential Coexpression

Differential coexpression analysis usually requires the definition of “distance” or “similarity” between measured datasets, the most common choice being Pearson’s correlation coefficients. However, Pearson’s correlation coefficient is sensitive to outliers [13]. Biweight midcorrelation is considered to be a good alternative to Pearson’s correlation since it is more robust to outliers. Example of a gene expression matrix is as follows:(1)Gene1Gene2⋯GenepZ1X11X12⋯X1pZ2X21X22⋯X2p⋮⋮⋮⋱⋮ZnXn1Xn2⋯Xnp.For each sampleZi, we measure expression levels for a set of genes, so Xij is the measurement of the expression level of the jth gene for the ith sample, where j=1,…,p. The xth column vector of matrix represents gene expression profile of gene X. In order to define the biweight midcorrelation (bicor) [13] of two numeric vectors x=(x1,…,xm) and y=(y1,…,ym), we first define ui, vi with i=1,…,m:(2)ui=xi-medx9madx,vi=yi-medy9mady,where med(x) is the median of vector x, mad(x) is the median absolute deviation of vector x, mad(x) is the median of new numeric vector in which each number is absolute difference between original vector value and med(x); this leads us to the definition of mad(x) and weight wi for xi, which are(3)madx=medxi-medX,wix=1-ui22I1-ui,where the indicator I(1-ui) takes 1 if 1-ui>0 and 0 otherwise. Thus, the weight wi(x) is close to 1 if xi is close to med(x), approaches 0 when xi differs from by nearly  9mad(x), and is 0 if xi differs from med(x) by more than 9mad(x). An analogous weight wi(x) can be defined for yi. Given the weights, we can define biweight midcorrelation of x and y as(4)bicorx,y=∑i=1mxi-medxwixyi-medywiy·∑j=1mxj-medxwjx2hhhih·∑k=1myk-medywky2-1.It should be noted that the equations of biweight midcorrelation do not involve an explicit identification of outliers, and all elements whose weightwi=0 can be considered outliers. The user can also set up the maximum allowed proportion of outliers using the argument “maxPOutliers”; the “maxPOutliers” is interpreted as the maximum proportion of low and high outliers separately. For the value of bicor from −1 to 1, −1 represents the maximum negative correlation and 1 represents the maximum positive correlation. Zero represents irrelevant correlation.


## 2.2. The “Differential Coexpression Threshold” Strategy

We used biweight midcorrelation to measure every pair of genes in the gene expression dataset and get a gene coexpression matrix. The gene coexpression matrix is a square and symmetric matrixP whose rows and columns correspond to the genes and whose element Pij denotes the coexpression relationship between genes i and j. In this paper, we use AGN which represents gene coexpression adjacency matrix in normal conditions and AGD which represents gene coexpression adjacency matrix in disease condition. To find differential coexpression gene modules which are coexpressed in normal condition and not related to disease condition, we set two thresholds T1 for adjacency matrix AGN in normal condition and T2 for adjacency matrix AGD in disease condition. AGN(i,j) is set to 1 if value of AGN(i,j) is greater than or equal to T1; otherwise, AGN(i,j) is set to 0 and AGD(i,j) is set to 1 if value of AGD(i,j) is less than or equal to T2; otherwise, AGD(i,j) is set to 0. We integrated AGN and AGD into a matrix AG after we had intersection of the corresponding elements of AGN and AGD. AG(i,j)=1 means coexpression value of gene i and gene j in AGN is greater than or equal to T1, and coexpression value of genes i and j in AGD is less than or equal to T2. AGN(i,j) also can be set to 1 if value of AGN(i,j) is less than or equal to T1 and AGD(i,j) is set to 1 if value of AGD(i,j) is greater than or equal to T2. The method is shown in (5).With the above mentioned strategy, we also setAG(i,j)=1 if the absolute value of AGN(i,j) subtracting AGD(i,j) is greater than or equal to T3 and the absolute value of AGN(i,j) is greater than or equal to the absolute value of AGD(i,j) simultaneously. This is a special type of coexpression change. In reality, coexpression reversal probably has biological significance. The coexpression reversal between normal condition and disease condition has advantage in disease. For example, the coexpression ofp53 andKlf4 recently reported that the positive or negative correlation between these two genes determines the outcome of DNA damage, DNA repair, or apoptosis [19]. We believe that our attention to this special coexpression change will help to explore subtle mechanisms involved in genes transcriptional regulation. We excavated maximum cliques which have biological significance from AG adjacency matrix to further investigate gene regulatory networks. Consider the following:(5)ifAGNi,j≥T1,thenAGNi,j=1,elseAGNi,j=0;ifAGDi,j≤T2,thenAGDi,j=1,elseAGDi,j=0;AGi,j=AGNi,j∩AGDi,jifAGNi,j≤T1,thenAGNi,j=1,elseAGNi,j=0;ifAGDi,j≥T2,thenAGDi,j=1,elseAGDi,j=0;AGi,j=AGNi,j∩AGDi,j.


## 2.3. The Threshold Selection Strategy

The two real value adjacency matrixes are transformed into a binary matrix which contains two elements 0 and 1 only. Choosing different thresholds will lead to different results; too largeT1 threshold or too small T2 threshold will lead to small link number, low density clique, and lost biological significance cliques. On the other hand, too small T1 or too large T2 will lead to many overlapping cliques. They are not helpful for finding biological significance differential coexpression gene disease-related modules. In fact, how to choose a reasonable threshold in conversion process is a problem which needs to be further studied. Generally, the selection of the threshold can be based on the proportion of outliers in the figure or the density of graph. The outlier is the point which is not connected to any edges. The density is defined as the ratio of number of edges to the maximum possible number of edges in the graph. The density of clique is 1.For gene expression data analysis, closely linked functional module is not the strict sense of maximum clique due to the lack of certain section. In this paper, we use density to measure approximation degree of functional module with gene differential coexpression clique, which may be having more biological significance.


## 2.4. The Maximum Clique Concept andk-Clique Algorithm

Graph theoretical concepts are useful for the description and analysis of interactions and relationships in biological systems. In gene coexpression graph, gene is represented by vertex and coexpression relationship by edge.G=(V,E) is an arbitrary undirected and weighted graph unless otherwise specified in graph theoretical concepts. V={1,2,…,n} is the vertex set of G, and E is the edge set of G. For each vertex i∈V, a positive weight wi is associated with i. AG=aijn×n is the adjacency matrix of G, where aij=1 if (i,j)∈E is an edge of G, and aij=0 if (i,j)∉E. Genes and relationship between genes are represented by vertex and edge, respectively.A graphG=(V,E) is complete if all its vertices are pairwise adjacent; that is, for all i,j∈V, (i,j)∈E. A clique C is a subset of V such that G(C) is complete. The maximum clique problem asks for a clique of maximum weight. An independent set (stable set and vertex packing) is a subset of V, whose elements are pairwise nonadjacent. The maximum independent set problem asks for an independent set of maximum cardinality. The size of a maximum independent set is the stability number of G (denoted by α(G)). The maximum weight independent set problem asks for an independent set of maximum weight. A maximum clique means a clique which is a subset of the nodes in V in which every pair of nodes in the subset is joined by an edge and is not a proper subset of any other cliques [20].In application, the identification of maximal cliques is often of limited interest since the requirement of complete connectivity is so restrictive. When dealing with imperfect systems or with experimental data, we may need to consider more general notions of cohesive subgroups. In this paper, we consider different notions of cohesive subgroups that includen-cliques, k-plexes, and λ-sets [18]. It is well known that the nodes of large real networks have a power law degree distribution [21]. Most real networks typically contain parts in which the nodes (units) are more highly connected to each other compared to the rest of the network. The sets of such nodes are usually called clusters, communities, cohesive groups, or modules [22–26], which have no widely accepted unique definition. The basic observation on which our modules definition relies is that a typical gene differential coexpression module consists of several complete (fully connected) subcliques that tend to share many of their nodes. To find meaningful communities, several basic requirements should be satisfied: it cannot be too restrictive, should be based on the density of links, is required to be local, should not yield any cut-node or cut link (whose removal would disjoin the community), and, of course, should allow overlaps. We employ the community definition specified above because none of the others in the literature satisfy all these requirements simultaneously [27–29].k-clique algorithm for detecting gene differential coexpression modules in a network has been published in the paper [26].k-clique algorithm is also named clique percolation method. The existing divisive and agglomerative methods recently used for large real networks have some disadvantages. Divisive methods cut the network into smaller and smaller pieces; each node is forced to remain in only one community and be separated from its other communities, most of which then necessarily fall apart and disappear [27, 30]. The agglomerative [31] method has the same problem. Thek-clique algorithm has demonstrated the advantages over the divisive method and agglomerative method. In the algorithm, although the numerical determination of the full set ofk-clique communities is a polynomial problem, the algorithm is exponential and significantly more efficient for the graphs corresponding to actual data. Thek-clique algorithm first locates all cliques (maximal complete subgraphs) of the network and then identifies the communities by carrying out a standard component analysis of the clique-clique overlap matrix [28]. Thek-clique algorithm uses the threshold probability d(k) (critical point) ofk-clique percolation to find all maximal complete subgraphs. The critical point is shown in (6), where N is the number of genes or vertex of graph:(6)dk=1k-1N1/k-1.Thek-clique algorithm gives two plausible choices to measure the size of the largestk-clique percolation cluster in (7) and (8). The most natural one, which we denote by N∗, is the number of vertices belonging to this cluster. ϕ is an order parameter associated with this choice as the relative size of that cluster:(7)ϕ=N∗N.The other choice is the numberL∗ ofk-cliques of the largestk-clique percolation cluster. The associated order parameter is again the relative size of this cluster:(8)φ=L∗L,where L denotes the total number ofk-cliques in the graph. L can be estimated as(9)L≈Nkdkk-1/2≈Nkk!dkk-1/2.In this paper, we use the biweight midcorrelation for constructing binary networks. Two-condition coexpression adjacency networks can always be transformed into a binary one by ignoring any directionality in the links and keeping only those stronger than a threshold weight. Then, the concept of maximum clique andk-clique algorithm were used to find gene differential coexpression modules. We named the proposed method “BMKC” (biweight midcorrelation andk-clique algorithm) method. Changing the threshold is like changing the resolution with which the community structure is investigated: by increasing, the communities start to shrink and fall apart. A very similar effect can be observed by changing the value of k as well: increasing k makes the communities smaller and more disintegrated but, at the same time, also more cohesive. More details aboutk-clique algorithm can be found in [28, 32].


## 3. Results

### 3.1. Experiment Result on Simulated Datasets

We first evaluate the algorithm in a supervised setting. We generate a control group of 30 samples and a disease group of another 30 samples, both consisting of 120 genes. For the control group, 20 coexpressed genes are sampled directly from the biweight midcorrelation. We focus on whetherk-clique algorithm can find coexpression gene modules from the background of noise. We first draw a vector with 20 rows and a vector with 30 columns from a standard normal distribution. The actual expression levels are obtained by adding independent errors sampled from a normal distribution with mean zero and standard deviation (SD) σ. These 20 genes form the target pattern. We then hide them in 100 additional noise genes, which are sampled independent and identically distributed (i.i.d.) from a standard normal distribution. The disease group is simulated by 120 independent noise genes drawn from a standard normal only.In the above setting, we use SDσ to tune the strength of the signal resulting from the 20 coexpressed genes. To observe its effect in detail, we use three different values: for a clear signal, σ=1/10, for medium noise, σ=1/4, and, for high noise, σ=1. To guard for sampling effects, we repeat each procedure 50 times and average the results, which are displayed in Figure 1. One can see that, for the clear and medium signal, the algorithm can recover the differentially coexpressed genes modules reliably. Also, depending on the prominence of the signal, the influence of σ is more or less pronounced. In an exploratory analysis setting with several hidden patterns, we could use T1, T2, and T3 to control the size of target patterns.Figure 1
The column bar graph shows the effect of the noise parameterσ on the size of the gene group found by our algorithm.


### 3.2. Analyzing a Type 2 Diabetes (T2D) in Rats

As a real-world application, we apply the BMKC method to a pair of type 2 diabetes (T2D) rats datasets (dataset pairT), which has been published in study [33]. Dataset pair T is from dataset GSE3068 of Gene Expression Omnibus (GEO) database. Yu et al. preprocessed dataset GSE3068. Dataset pair T includes 4765 genes in 10 disease samples and 10 normal samples. We use our algorithm to find differential coexpression modules in the type 2 diabetes.For computational efficiency, we calculate the sum of each row or column of adjacency matrix; the sum means the number of genes related to the gene. The gene is outlier if the sum is zero. First, we calculate the sum of each row or column of the adjacency matrix and delete the outlier. Second, we calculate the sum of each row or column of the adjacency matrix and discard the lower 50% of them. We setT3=1.3 and the minimum number of each clique to four. Finally, we apply our algorithm to the remainder genes and excavate two differential coexpression modules. Tables 1 and 2 list each gene symbol in the clique. The adjacency graphs of each differential coexpression module are shown in Figures 2 and 3. From these two figures, we can see that the cliques in each of the differential coexpression modules are overlapping, forming a closely related module. In normal condition, the absolute bicor value of total of 24 genes in modules distributes from 0.78 to 0.97. Yet, in disease condition, the absolute bicor value of genes distributes from 0.21 to 0.09. In the results of our study, the gene differential coexpression modules included quite a number of previously reported T2D-related genes:Hifla andSirt2 [34],Smarca4 [35],Sh2b2 [36],Madd [37], andRxrb [38]. Despite not being previously reported to be related with T2D, other genes in the modules should receive adequate attention for their distinct traits from the perspective of differential coexpression. Further studies on the transcriptional mechanisms and functional consequences could pay more attention to these genes.Table 1
Genes in each clique.
Clique number
Gene symbol
1
Hifla
*
Ifngr1
RGD1305094
Tenc1
Sirt2
2
Clcn1
Smarca4
Zkscan17
Rpl27a
Sirt2
3
Hifla
Ifngr1
Pfkfb3
Tenc1
Sirt2
4
Sh2b2
Pcsk5
Lamc1
Rpl27a
Sirt2
5
Lamc1
Smarca4
Zkscan17
Sirt2
Rpl27a
6
Hifla
RGD130504
Mxd4
Sirt2
7
Tra1
Smarca4
Zkscan17
Sirt2
∗Bold genes refer to the previously reported T2D-related genes. The other genes are identified in the differential coexpression modules.Table 2
Genes in each clique.
Clique number
Gene symbol
1
Vars2
Apba3
Madd
Zkscan17
Ppil3
Rxrb
2
Snd1
Madd
Pctk2
RxrbFigure 2
The adjacency graph of first gene differential coexpression module.Figure 3
The adjacency graph of second gene differential coexpression module.


### 3.3. Significance Analysis of the BMKC Method

Naturally, the question of whether our findings are artifacts of the high dimensionality of the data arises. To assess this question, we apply a permutation procedure. Under the null hypothesis, we assume that all genes are mutually independent in both conditions groups. We heuristically sample from the null hypothesis by (group-wise) shuffling the expression values for each gene independently. Thus, random expression data are generated where all covariance structures are removed. Applying our algorithm to the randomized data yields one random score. We repeat the procedure 1000 times. Using the empirical distribution of the simulated scores, the simulated score means the global total sum of differential coexpression change of each gene in modules. We calculateP values for the observed scores in the nonpermuted data. For each of the patterns in the type 2 diabetes example, we only observe one random score smaller than the biological one. This corresponds to an empirical P value of 0.001. Hence, it is unlikely that the observed differential coexpression is a chance artifact.


## 3.1. Experiment Result on Simulated Datasets

We first evaluate the algorithm in a supervised setting. We generate a control group of 30 samples and a disease group of another 30 samples, both consisting of 120 genes. For the control group, 20 coexpressed genes are sampled directly from the biweight midcorrelation. We focus on whetherk-clique algorithm can find coexpression gene modules from the background of noise. We first draw a vector with 20 rows and a vector with 30 columns from a standard normal distribution. The actual expression levels are obtained by adding independent errors sampled from a normal distribution with mean zero and standard deviation (SD) σ. These 20 genes form the target pattern. We then hide them in 100 additional noise genes, which are sampled independent and identically distributed (i.i.d.) from a standard normal distribution. The disease group is simulated by 120 independent noise genes drawn from a standard normal only.In the above setting, we use SDσ to tune the strength of the signal resulting from the 20 coexpressed genes. To observe its effect in detail, we use three different values: for a clear signal, σ=1/10, for medium noise, σ=1/4, and, for high noise, σ=1. To guard for sampling effects, we repeat each procedure 50 times and average the results, which are displayed in Figure 1. One can see that, for the clear and medium signal, the algorithm can recover the differentially coexpressed genes modules reliably. Also, depending on the prominence of the signal, the influence of σ is more or less pronounced. In an exploratory analysis setting with several hidden patterns, we could use T1, T2, and T3 to control the size of target patterns.Figure 1
The column bar graph shows the effect of the noise parameterσ on the size of the gene group found by our algorithm.


## 3.2. Analyzing a Type 2 Diabetes (T2D) in Rats

As a real-world application, we apply the BMKC method to a pair of type 2 diabetes (T2D) rats datasets (dataset pairT), which has been published in study [33]. Dataset pair T is from dataset GSE3068 of Gene Expression Omnibus (GEO) database. Yu et al. preprocessed dataset GSE3068. Dataset pair T includes 4765 genes in 10 disease samples and 10 normal samples. We use our algorithm to find differential coexpression modules in the type 2 diabetes.For computational efficiency, we calculate the sum of each row or column of adjacency matrix; the sum means the number of genes related to the gene. The gene is outlier if the sum is zero. First, we calculate the sum of each row or column of the adjacency matrix and delete the outlier. Second, we calculate the sum of each row or column of the adjacency matrix and discard the lower 50% of them. We setT3=1.3 and the minimum number of each clique to four. Finally, we apply our algorithm to the remainder genes and excavate two differential coexpression modules. Tables 1 and 2 list each gene symbol in the clique. The adjacency graphs of each differential coexpression module are shown in Figures 2 and 3. From these two figures, we can see that the cliques in each of the differential coexpression modules are overlapping, forming a closely related module. In normal condition, the absolute bicor value of total of 24 genes in modules distributes from 0.78 to 0.97. Yet, in disease condition, the absolute bicor value of genes distributes from 0.21 to 0.09. In the results of our study, the gene differential coexpression modules included quite a number of previously reported T2D-related genes:Hifla andSirt2 [34],Smarca4 [35],Sh2b2 [36],Madd [37], andRxrb [38]. Despite not being previously reported to be related with T2D, other genes in the modules should receive adequate attention for their distinct traits from the perspective of differential coexpression. Further studies on the transcriptional mechanisms and functional consequences could pay more attention to these genes.Table 1
Genes in each clique.
Clique number
Gene symbol
1
Hifla
*
Ifngr1
RGD1305094
Tenc1
Sirt2
2
Clcn1
Smarca4
Zkscan17
Rpl27a
Sirt2
3
Hifla
Ifngr1
Pfkfb3
Tenc1
Sirt2
4
Sh2b2
Pcsk5
Lamc1
Rpl27a
Sirt2
5
Lamc1
Smarca4
Zkscan17
Sirt2
Rpl27a
6
Hifla
RGD130504
Mxd4
Sirt2
7
Tra1
Smarca4
Zkscan17
Sirt2
∗Bold genes refer to the previously reported T2D-related genes. The other genes are identified in the differential coexpression modules.Table 2
Genes in each clique.
Clique number
Gene symbol
1
Vars2
Apba3
Madd
Zkscan17
Ppil3
Rxrb
2
Snd1
Madd
Pctk2
RxrbFigure 2
The adjacency graph of first gene differential coexpression module.Figure 3
The adjacency graph of second gene differential coexpression module.


## 3.3. Significance Analysis of the BMKC Method

Naturally, the question of whether our findings are artifacts of the high dimensionality of the data arises. To assess this question, we apply a permutation procedure. Under the null hypothesis, we assume that all genes are mutually independent in both conditions groups. We heuristically sample from the null hypothesis by (group-wise) shuffling the expression values for each gene independently. Thus, random expression data are generated where all covariance structures are removed. Applying our algorithm to the randomized data yields one random score. We repeat the procedure 1000 times. Using the empirical distribution of the simulated scores, the simulated score means the global total sum of differential coexpression change of each gene in modules. We calculateP values for the observed scores in the nonpermuted data. For each of the patterns in the type 2 diabetes example, we only observe one random score smaller than the biological one. This corresponds to an empirical P value of 0.001. Hence, it is unlikely that the observed differential coexpression is a chance artifact.


## 4. Conclusions

In this paper, we proposed a new approach in gene sets level for differential coexpression analysis, which combine biweight midcorrelation and threshold selection strategy and also applied maximum clique concept withk-clique algorithm to the specific gene set to further investigate gene regulatory networks. Biweight midcorrelation is more robust for outliers and threshold selection strategy is an effective preprocess step of the proposed method. Experimental results on simulated datasets show that our method had good performance. We apply the proposed BMHT method to real dataset designed for T2D study, and two differential coexpression gene modules were detected, which should be a useful resource for T2D study and could be used for exploring the biological function of the related genes. In the future, we will focus on how to quickly excavate gene differential coexpression module from gene coexpression adjacency matrix.



---

*Source: 836929-2015-08-03.xml*