 Methodology Article
 Open access
 Published:
Assisted clustering of gene expression data using ANCut
BMC Genomics volume 18, Article number: 623 (2017)
Abstract
Background
In biomedical research, gene expression profiling studies have been extensively conducted. The analysis of gene expression data has led to a deeper understanding of human genetics as well as practically useful models. Clustering analysis has been a critical component of gene expression data analysis and can reveal the (previously unknown) interconnections among genes. With the high dimensionality of gene expression data, many of the existing clustering methods and results are not as satisfactory. Intuitively, this is caused by “a lack of information”. In recent profiling studies, a prominent trend is to collect data on gene expressions as well as their regulators (copy number alteration, microRNA, methylation, etc.) on the same subjects, making it possible to borrow information from other types of omics measurements in gene expression analysis.
Methods
In this study, an ANCut approach is developed, which is built on the regularized estimation and NCut techniques. An effective R code that implements this approach is developed.
Results
Simulation shows that the proposed approach outperforms direct competitors. The analysis of TCGA (The Cancer Genome Atlas) data further demonstrates its satisfactory performance.
Conclusions
We propose a more effective clustering analysis of gene expression data, with the assistance of information from regulators. It provides a new venue for analyzing gene expression data based on the assisted analysis strategy.
Background
In biomedical research, profiling studies have been extensively conducted. Information so collected has led to a better understanding of human genetics as well as practically useful models. In genetic research, gene expression data have been playing an essential role in the past decades. Compared to DNA and epigenetic changes, gene expressions are “closer” to phenotypes. With relatively mature techniques, they are also easy to measure at a genomewide scale. Extensive methodological studies have been conducted on how to more effectively analyze gene expression data.
In the analysis of gene expression data, clustering has been playing an essential role. In some studies, clustering has been used to suggest/identify unknown functions of genes, with genes in the same cluster likely having related biological functions [1, 2]. In some other studies, clustering has been used as a way of reducing dimensionality. For example, some studies have suggested that conducting principal component analysis (and other analysis) on genes within the same clusters to reduce dimensionality is more sensible than doing that on all genes [3, 4]. In the literature, a large number of clustering methods have been developed and applied to gene expression data. Examples include the Kmeans, hierarchical clustering, agglomerative clustering, graphbased clustering, modelbased sparse clustering, and others. For reviews and comprehensive discussions, we refer to [5, 6]. Although having certain technical differences, most of these methods share the common strategy of reducing “distance” within clusters while maximizing “distance” across clusters.
Despite great successes, it has also been suggested that clustering analysis of gene expression data still very often generates unsatisfactory results [7]. Gene expression studies usually have a limited sample size but a moderate to large number of genes. Intuitively, there is a lack of “sufficient information”. Multiple strategies have been developed to increase information. One strategy is to increase sample size via pooled analysis [8], which demands additional samples. There are also studies that use additional, especially biological, information [9]. In this study, we develop a different strategy which takes advantage of additional omics measurements on the same samples. Gene expression levels are regulated by multiple mechanisms, including copy number alteration, microRNA, methylation, and others. Intuitively, such regulators contain additional information on gene expressions. In recent biomedical research, a prominent trend is to conduct “multidimensional” profiling and collect data on gene expressions as well as their regulators on the same samples, making it possible to “borrow information” from other types of omics measurements for gene expression data analysis.
In this study, our goal is to conduct more effective clustering analysis of gene expression data. With the assistance of information from regulators. A novel clustering approach is developed to achieve this goal. The basic strategy is similar to that of minimizing distance within clusters and maximizing distance between clusters, and so the proposed approach has a solid statistical ground. Advancing from the existing methods, information of regulators is utilized in clustering gene expression data. More accurate clustering can be achieved by using more information. This is especially desirable with the decreasing cost of profiling but increasing cost of sample collection. This study differs from most of the existing ones that conduct the integrated analysis of multidimensional omics data. Studies such as [10, 11] aim at building more accurate disease outcome models by integrating gene expression and regulator data. Studies such as [12] aim at identifying “hot spots” of the chromosome that host multiple types of omics changes. In contrast, our goal is to conduct clustering, which is a more fundamental goal of gene expression data analysis. This study is warranted with an important analysis goal and an innovative and effective new method.
In what follows, we first describe the data structure. For the simplicity of notation, we use copy number alteration (CNA) as a representative of gene expression (GE) regulators. The proposed method will be directly applicable to other types of regulators. We then describe the proposed method, computational algorithm, and software development. Simulation and the analysis of TCGA (The Cancer Genome Atlas) data are conducted to gauge performance of the proposed approach relative to the alternatives. It is noted that although developed for GE data, the proposed approach can have broader applications.
Methods
Consider a dataset with n iid samples. For the ith sample, assume that measurements are available on p GEs, denoted as Y _{ i }=(Y _{ i1},Y _{ i2},⋯,Y _{ ip })^{′}. In addition, assume that measurements are also available on q CNAs, denoted as X _{ i }=(X _{ i1},X _{ i2},⋯,X _{ iq })^{′}. For another type of regulator (for example, microRNA), the proposed approach is directly applicable. If there are two or more types of regulators, following [13], we can stack them together and create a “mega” vector of regulators. Our strategy is to use information in X _{ i }’s to assist the clustering of Y _{ i }’s. In the next subsection, we first describe our strategy for modeling the regulation between GEs and regulators. The assisted clustering method will then be developed.
Modeling the GECNA regulation
Following the literature [14, 15], we describe the GECNA regulation using regression. Specifically, consider
where β is the matrix of unknown regression coefficients, and ε _{ i } is the vector of “random errors”. In the literature, there are multiple ways of modeling the GECNA regulation. The regression approach has been adopted in quite a few recent studies and shown to have advantages over many alternatives for example the correlationbased. It is especially suitable for analysis with a large number of GEs and CNAs.
The regulation relationship is reflected in β, with a nonzero component corresponding to a regulation between a GE and a CNA and the magnitude describing the strength of regulation. For estimating β, we consider the penalized estimate
where Y and X are matrices consisting of Y _{ i }’s and X _{ i }’s, and λ>0 and 0≤α≤1 are datadependent tuning parameters. The penalization approach is adopted to accommodate the high data dimensionality and for selection: for a specific gene, its expression level is expected to be affected by only a few CNAs, and a CNA is expected to affect the expression levels of only a few GEs, which poses a variable selection problem. The elastic net (Enet) penalty is adopted for its simplicity and to accommodate (possibly high) correlations among CNAs [16]. Similar estimation approaches have been adopted in the literature [17]. In data analysis, this estimation is effectively realized using the R package glmnet. The two tuning parameters λ and α are selected using Vfold cross validation (V=5 in our numerical study).
With the estimate \(\hat {\boldsymbol {\beta }}\), denote the “predicted” GE values as \(\hat {\boldsymbol {Y}}=\hat {\boldsymbol {\beta }} \boldsymbol {X}\), which describe the component of GEs regulated by the regulators. Accordingly, \(\boldsymbol {Y}^{c}= \boldsymbol {Y} \setminus \hat {\boldsymbol {Y}}\) contains the levels of GEs regulated by other regulators (that are not included in X), affected by other mechanisms, as well as “random variations”. This decomposition strategy for GEs has been recently developed in the literature [18] under other contexts and shown to provide important additional insights into GEs (beyond treating GEs as a whole). More discussions are also provided in the next section.
Assisted clustering
For GEs, consider the weight matrix W=(w _{ jl })_{ p×p }, where the nonnegative element w _{ jl } measures the similarity between genes j and l. For a pair of the original GE measurements (as included in Y), we define w _{ jl } equal to the inverse of their Euclidean distance. Note that there are multiple ways of defining the similarity. This definition is adopted because of its simplicity. It shares a similar spirit with the popular Kmeans approach. Further, we define \(\widehat {\boldsymbol {W}}\), which is obtained in a similar way as W but using \(\hat {\boldsymbol {Y}}\), the regulated component of GEs.
Denote A _{1},…,A _{ K } as a partition of {1,…,p} which leads to K disjoint clusters. For A _{ k }, denote \(A_{k}^{c}\) as its complement. We propose the ANCut (Assisted NCut) measure as
where
and
With a fixed K, the optimal clustering minimizes the ANCut measure.
Rationale The proposed approach has been motivated by the following considerations. It is built on the NCut technique, which is originally developed in imaging and other scientific fields [19] and more recently applied to genetic and other data types [20, 21]. The NCut technique may have multiple advantages over the alternatives. Specifically, the “cutting” step is relatively independent of the similarity/distance construction. Without making restrictive assumptions on the similarity measure and underlying data distributions and models, it enjoys very broad applicability. In addition, both the numerator and denominator have lucid interpretations, with the numerator measuring the acrosscluster similarity and the denominator measuring the withincluster similarity. More discussions on the NCut technique are provided in Additional file 1.
The most significant advancement from the existing approaches is that the numerator and denominator in (2) are defined using two different sets of similarity measures. The levels of GEs can be affected by multiple mechanisms. (a) One is regulating by the regulators measured in X. (b) Many (or most) existing studies are not “exhaustive” and do not measure all regulators. Thus, possibly there are regulators not included in X. (c) There are mechanisms other than regulation that may also affect GE levels. For example, the expression level of a gene can be affected by other genes, for example through RNA interference (which is also known as cosuppression or posttranscriptional gene silencing). Also, one gene can code for a transcription factor, and it can bind to the promoter region of another gene, which consequently affects its expression level. Following the literature [13, 18], in our analysis, we decompose GE levels into two components: the first is (a), and the second is (b)+(c).
Published studies [22, 23] under other contexts have shown that jointly considering the two components of GEs can lead to more sensible analysis results than considering them as a whole. Motivated by such observations, our strategy is to minimize the NCut measures for both components of GEs. Since the overall GE levels are equal to the sum of the two components, loosely speaking, using W (which measures the sum) and \(\widehat {\boldsymbol {W}}\) is equivalent to using the two individual components. Our numerical analysis also confirms this intuition (results omitted). With this strategy, a seemingly more “straightforward” and “symmetric” objective function is \(\sum \limits _{k=1}^{K}\frac {\text {cut}\left (A_{k},A_{k}^{c};\boldsymbol {W}\right)} {\text {cutvol}\left (A_{k};\widehat {\boldsymbol {W}}\right)}+\sum \limits _{k=1}^{K}\frac {\text {cut}\left (A_{k},A_{k}^{c};\widehat {\boldsymbol {W}}\right)} {\text {cutvol}\left (A_{k}; \boldsymbol {W}\right)}.\) Our exploration suggests that this objective function is computationally more expensive but leads to similar results as the one in (2).
A toy example To demonstrate the operating characteristics of the proposed method, we consider a toy example with 10 GEs and 10 CNAs. Data generation is the same as described in detail in the next section, except with a lower dimensionality. The true data structure is shown in the left panel of Fig. 1. There are a total of two equalsized GE clusters, represented by different colors. The degree of similarity between two GEs (and CNAs) is represented by the thickness of lines. For the simulated dataset presented in Fig. 1, the proposed analysis is able to fully recover the true data structure with the penalized estimation and the true clustering structure with the ANCut approach. As an alternative, we also consider the popular Kmeans approach. We observe that the Kmeans approach completely fails for this dataset. It identifies two clusters one with nine GEs and the other with one GE. We have also experimented with some other alternatives and observe similar failing performance.
Computation
For optimizing the objective function defined in (2), we adopt the simulated annealing (SA) technique [24]. Use t as the index of iteration. At iteration t, denote \(A^{(t)}=\left \{A_{1}^{(t)},\dots,A_{K}^{(t)}\right \}\) as the partition (clustering result) and ANCut(t) as the value of the objective function. Further denote B as the maximum number of iterations. The value of B is not important, as long as it is large enough. Define the temperature function as T(t)=Llog(t+1). In our numerical study, we set L=1000, which generates satisfactory results. In practice, to be prudent, other L values may also need to be examined. Extensive discussions on tuning parameter selection with the SA technique are available in the literature. The proposed algorithm proceeds as follows.
Step 1 Randomly initialize \(A^{(0)}=\left \{A_{1}^{(0)},\ldots,A_{K}^{(0)}\right \}\). In our numerical study, different initial values lead to almost identical results.
Step 2 Set t=t+1. For k=1,…,K, compute p _{ i } as the number of (j,l) pairs with \(j, l \in A_{k}^{(t1)}\). Draw k(−) and k(+) from {1,…,K} with probabilities proportional and inversely proportional to p _{ i }, respectively.
Step 3 Draw i randomly from \(A_{k()}^{(t)}\). Set \(A_{k(+)}^{(t)}=A_{k(+)}^{(t1)}\cup \{i\}\) and \(A_{k()}^{(t)}=A_{k()}^{(t1)}\setminus \{i\}\). For j≠k(+),k(−), \(A_{j}^{(t)}:=A_{j}^{(t1)}\).
Step 4 If ANCut(t)≤ANCut(t−1), keep A ^{(t)} as it is. If not, keep A ^{(t)} as it is with probability \(\exp \Bigg (\frac {\text {ANCut} (t) \text {ANCut} (t1)}{T(t)} \Bigg)\), and otherwise A ^{(t)}=A ^{(t−1)}.
Step 5 Repeat Steps 24 until t=B.
Extensive research on the SA technique is available in the literature [24]. In Step 2, the proposed probabilities prefer adding a new member to a small cluster and deleting a member from a large cluster. Thus, the “prior” is that clusters have similar sizes. Note that this strategy may be somewhat subjective and can be adjusted according to the specific scientific context. Convergence of the SA algorithm to the global optimizer has been established in the literature [24]. It is achieved in all of our numerical examples.
The proposed analysis, which consists of penalized estimation and ANCut, is computationally affordable. The two steps have computational complexity O(n p q) and O(B p), respectively. For a simulated dataset with p=q=200 and n=50, we consider 100 tuning parameter values in penalized estimation and B=10,000 in ANCut. The proposed analysis takes about 30 s on a laptop with standard configurations.
Software To facilitate data analysis, we develop an R package NCutYX and make it publicly available at https://github.com/shuanggema. If the R library devtools is installed, then the NCutYX package can be easily installed using devtools::install_github(“shuanggema/NCutYX”). The proposed approach is implemented using the function ANCut, which proceeds as follows: clust ←ANCut(Y, X, K = 2, B = 3000, L = 1000, alpha = 0.5, nlambdas = 100, ncv = 5, dist = “euclidean”) In the above command, Y is the data matrix of GEs, X is the data matrix of regulators, K is the number of clusters, B is the number of SA iterations, L is the temperature coefficient, alpha is α in the Enet penalty (note that in this sample command, we fix α as its default value. It is easy to datadependently select α), nlambdas is the number of λ values in Enet, ncv is the number of crossvalidations, and dist specifies that the Euclidean distance is used in defining the dissimilarity. The resulting object clust is a list where the first entry (clust[[1]]) is a vector of SA sequence, the second entry (clust[[2]]) includes the clustering results, and the third entry (clust[[3]]) contains the optimal λ value. Written in a friendly way, the package can be easily adopted and modified.
Results
Simulation
Use the same notations as in the last section. In simulation, GEs and CNAs are linked with the regression model Y=β X+ε. Each element of the matrix ε, denoted by ε _{ i,j }, is iid from a normal distribution with mean=0 and sd=2. For each subject, the CNA values are generated from a multivariate normal distribution with marginal means 0 and variancecovariance matrix Σ. Σ has a blockdiagonal structure with two blocks, each of which has size q/2. In each block, the diagonal elements are equal to 1, and all offdiagonal elements are equal to ρ. In TCGA (which is analyzed in the next section) and other datasets, it has been observed that the processed CNA data have unimodal distributions close to normal (although the raw data may have different distributional characteristics). Two equalsized CNA clusters are generated, with those in the same cluster correlated and different clusters uncorrelated. For the regression coefficient matrix, consider
For each column, with the sparsity consideration, we set q _{0} coefficients to be nonzero and the rest zero. We consider two sparsity levels with q _{0}= 3 and 6. Two different coefficient settings are considered. The first (C1) has the nonzero coefficients randomly generated from Uniform [h/2,h], where the parameter h determines the strength of regulation. The second (C2) has the nonzero coefficients randomly generated from Uniform [−0.15,0.25], which describes the scenario where CNAs have both positive and negative regulations on GEs and the amount of positive and negative regulation is not equal. Under this data generating structure, there are two clusters of CNAs/GEs, and GEs in the first (second) cluster are regulated by CNAs in the first (second) cluster. Beyond sparsity, this setting also describes the “localization” of regulations.
Under this simulation setting, the component of GE that is regulated by regulators, i.e., effect “(a)” in the “Assisted clustering” section, is β X. For the other component of GE, i.e., effects “(b)+(c)”, is randomly generated. There are at least two considerations for this. First, unlike for the regulators’ effects, research on the second component of GE is still much limited. It is not entirely clear how to simulate effects “(b)+(c)”. More importantly, this setting favors the alternative approaches (to be described below). If the proposed approach has competitive performance under this unfavorable setting, it is reasonable to expect better results under favorable settings.
When evaluating performance of the proposed approach, we consider both accuracy and stability. With a set of clusters {A _{1},…,A _{ K }}, an adjacency matrix A=(a _{ jl })_{ p×p } can be constructed, where the (j,l)th element a _{ jl }=1 if j,l∈A _{ k },k=1,⋯,K and 0 otherwise. Let A _{ T } and \(\hat {\boldsymbol {A}}\) be the adjacency matrices of the true and estimated clusters, respectively. Then the accuracy measure is defined as the diversity between A _{ T } and \(\hat {\boldsymbol {A}}\),
where ⊙ is the componentwise product. With N replicates, denote \(\widehat {\boldsymbol {A}}^{(1)},\ldots,\widehat {\boldsymbol {A}}^{(N)}\) as the adjacency matrices as constructed above. We define the stability measure as
This definition shares a similar spirit with the Ustatistic stability measure [25]. For M _{ accuracy }, a smaller value suggests a higher accuracy; For M _{ stability }, a smaller value suggests more stable results.
A closer look at the proposed approach
Consider the setting with coefficient C1, with p=100, and q=p/2. In addition, set ρ=0.10 and h=0.10. We first simulate one replicate and present the heatmaps for the observed and predicted GE levels in Fig. 2. It is observed that the heatmap using the predictive GEs shows a clearer clustering structure. This provides an intuitive justification for the proposed strategy of using CNA information.
Consider the true clustering structure (A _{1},A _{2}). Further consider \(\text {NCut}(A_{1}, A_{2})=\sum \limits _{k=1}^{2}\frac {\text {cut}\left (A_{k},A_{k}^{c};\boldsymbol {W}\right)} {\text {cutvol}\left (A_{k};\boldsymbol {W}\right)}\), which is the objective function defined in (2) using the observed GEs for both numerator and denominator, and ANCut(A _{1},A _{2}). With 100 simulated replicates, the mean values are
This suggests that the proposed assisted analysis can more strongly define the clustering structure.
We further examine the potential gain by using information of CNAs for clustering GEs. Specifically, beyond the proposed ANCut, we also consider the popular Kmeans approach. As the simulated GEs have normal distributions, Kmeans is favored. In addition, we also consider the approach ANCut_{ T }, which is the proposed ANCut but with β in the place of \(\hat {\boldsymbol {\beta }}\). Under this approach, we have perfect information on the GECNA regulation. As shown in Table 1, multiple scenarios on n, p, and h are considered. With 100 replicates, we compute the mean accuracy measure M _{ accuracy }. It is observed that if the regulation is perfectly known, then using ANCut leads to perfect identification of clusters. If the regulation needs to be estimated, using the proposed ANCut leads to accuracy better than or similar to that of Kmeans. This further justifies the value of using regulator information.
Comparison with the alternative methods
To better gauge performance of the proposed method, we compare with competing alternatives including Kmeans and spectral clustering (Spec.) [26]. These two alternatives are considered because of their popularity and satisfactory performance observed in published studies [27–29]. We also consider NCut to better appreciate the merit of assisted analysis. We comprehensively consider multiple combinations of n,p,q,q _{0},h, and ρ values for coefficient setting C1. We also examine multiple scenarios under coefficient setting C2 with various n,p,q,q _{0}, and ρ=0.40. There are a total of 80 simulation settings. Under each setting, we simulate 100 replicates. Summary statistics are presented in Tables 2, 3, 4, 5, and 6. Performance of the proposed approach depends on the strength of GECNA regulation, correlation of CNAs, data dimensionality, sample size, and others, in a similar way as observed for other clustering methods. Across the whole range of simulation settings, the proposed method is observed to have competitive performance, with superior accuracy and stability. For example, in Table 2 with n=200, p=500, q=500, and q _{0}=6, the proposed method has M _{ accuracy } equal to 28.4%, compared to 42.4% (NCut), 41.8% (Kmeans), and 37.4% (Spec.). It also has satisfactory performance in terms of stability. Under this specific setting, the M _{ stability } values are 17.6% (proposed), 40.6% (NCut), 42.2% (Kmeans) and 34.5% (Spec.). It is observed that performance of the proposed method decays with the increase of dimensionality and correlation among CNAs. This observation is reasonable as the proposed method involves estimating the regulation relationship. The matrix β has a total of q×p parameters, which can be very difficult to estimate with a moderate sample size. This estimation gets challenged with an increase in data dimensionality and correlation. When CNAs have both positive and negative regulation effects on GEs, of which the results are provided in Table 6, the proposed method also performs better than the alternatives. For example, with n=400, p=500, q=250, and q _{0}=6, the proposed method has M _{ accuracy } equal to 33.6%, compared to 34.5% (NCut), 39.4% (Kmeans) and 40.4% (Spec.). It also has satisfactory performance in terms of stability. Under this specific setting, the M _{ stability } values are 15.2% (proposed), 19.3% (NCut), 17.3% (Kmeans) and 16.2% (Spec.).
When data are available on both GEs and regulators, studies under other contexts [11, 30] have directly conducted their joint analysis and observed improvement (over analyzing GE only). Here we also briefly experiment with such an approach. Specifically, in Table 3, we also consider an approach that we call augmented Kmeans (AugK). This approach is based on the Kmeans and includes both GEs and CNAs in the clustering algorithm. Specifically, Kmeans is adopted to analyze the stacked data (Y,X). To compare with other approaches, only the cluster memberships of GEs are tracked in the final results. Simulation results in Table 3 suggests that, unlike in regression analysis [17], directly integrating GEs and regulators fails in clustering analysis. Specifically, it tends to cluster different data types together, as opposed to clustering connected GEs and CNAs together. It is noted that this approach seems very stable. This is also reasonable. It tends to cluster all GEs in one big cluster. Thus, even though the results are wrong, they are very stable.
Remarks
Beyond those described above, we have also experimented with a few other settings and observed similar satisfactory performance of the proposed approach. In the above simulation, we have set K=2. We have examined similar simulation settings with K=3,…,10 and made similar observations. We have compared with the most popular alternatives, whose stable and satisfactory performance has been well observed in the literature. It is of interest to conduct more extensive comparisons with other (less popular) approaches in future studies.
Data analysis
TCGA (https://tcgadata.nci.nih.gov/docs/publications/tcga/?) is one of the most recent multidimensional studies and has collected multiple types of omics measurements on the same subjects for multiple cancer types. The data have been recently collected and published and have a high quality. We analyze data on cutaneous melanoma (SKCM), which poses a serious public health concern. In this study, we analyze the processed level 3 data, which are downloaded from TCGA Provisional using the R package CGDS. We focus on metastatic samples of the Whites. Detailed information on data collection and processing is available on the TCGA website and elsewhere. Briefly, GE data were collected using the IlluminaHiseq RNAseq V2 platform and have been lowessnormalized, logtransformed, and mediancentered. The robust Zscores represent the gene expression status (up or down regulated) in tumor samples relative to normal tissues. A total of 19,626 measurements are available on 371 samples. CNA measurements were obtained using the Affymetrix Genomewide Human SNP array 6.0 platform. The loss and gain levels of copy number changes of tumors compared to normal tissues are identified using segmentation analysis and expressed in the log2 transformed form. A total of 21,699 measurements are available on 366 samples. GE and CNA data are merged using sample ID, resulting in a total of 366 samples. There are a few missing measurements, which are filled in using imputation. Performance of the proposed approach (and alternatives) decays when the ratio of the number of unknown parameters and sample size increases. We first conduct a simple prescreening via marginal analysis based on overall survival (an important clinical outcome) and select the top 1000 most significant genes. We then use GOTerm Finder [31, 32] and search for the GO (Gene Ontology) biological processes. 382 out of the 1000 genes have well defined GO terms and are selected for downstream analysis.
Unlike in simulation, the number of clusters is unknown. Determining the optimal number of clusters is a challenging problem. In our data analysis, we adopt the GAP approach [33], which has been coupled with multiple clustering methods and extensively adopted in data analysis. With the GAP approach, K=4 clusters are constructed with the 382 GEs. For Kmeans, for comparability, we also set K=4. The spectral clustering is not considered because of its inferior performance in simulation. With ANCut, the cluster sizes are 110, 105, 97, and 70, respectively. With Kmeans, the cluster sizes are 27, 350, 1, and 4, respectively. The Kmeans results are less desirable with the dominating majority of genes in a single cluster. The other alternatives also have similar unsatisfactory results. Detailed clustering results are available from the authors.
The 382 genes have a total of 104 GO biological processes. More closely examining these processes suggests that the majority are related to “regulation”. We further separate the 104 processes into four categories: positive regulation, negative regulation, regulation (without a well defined “direction”), and other. In Fig. 3, we compute the proportions of genes in the four clusters that have the four categories of processes. For ANCut, differences across the four clusters are clearly seen. For example, cluster 4 has a higher percentage of “regulation”, whereas cluster 1 has a higher percentage of “other”. With the highly unbalanced clustering, the Kmeans results, also presented in Fig. 3, are less interesting. For all of the 104 processes, we present their distributions in the four clusters in Fig. 4. Each bar represents one process, and different colors in the bar denote different proportions of genes in the clusters. Due to space limitation, the processes’ names are not shown. In general, the distributions of processes are quite different across clusters, which suggests the effectiveness of the proposed clustering. We further examine the ten most representative processes in detail. The results are presented in Fig. 5. It is shown that, for example, cluster 4 has a much higher percentage of “immune system process” than the other three clusters, while cluster 3 has a higher percentage of “localization in cell” than the others. Further examining the genes suggests that those annotated to the same processes are highly likely to be clustered together. For example, 90%, 88%, and 88% of the genes annotated to the processes “response to type I interferon”, “type I interferon signaling pathway”, and “cellular response to type I interferon” belong to cluster 4. With Kmeans, as the majority of genes are in a single cluster, such analysis is not conducted. The sensible biological findings provide support to the validity of the proposed clustering. Another finding from data analysis is that the results (distributions of biological processes) do depend on the number of clusters. The results with K=3 are presented in Figure A1 (Additional file 1). This finding is reasonable and has also been observed with other clustering methods in the literature.
Discussion
Clustering is an important step of gene expression data analysis. Beyond having independent value, it has also served as the foundation of many other analyses. Although a large number of methods have been developed in the literature, since the practical results are often unsatisfactory, there is still a great demand for more effective methods. In the analysis of genetic data, a major problem is the “lack of information”. To tackle this problem, we have proposed conducting assisted analysis by borrowing information from the regulators of GEs. The proposed method is built on the NCut technique, which has multiple notable advantages but still limited applications in genetic data analysis. The proposed ANCut has been partly motivated by the recent “decomposing gene expressions” strategy and is the first to do so in clustering analysis. It has an intuitive formulation and can be effectively realized using an SA algorithm. In simulation, it outperforms the direct competitors with better accuracy and stability. We acknowledge that there are many other possible simulation settings and potentially applicable alternatives. The considered simulation settings have comprehensively covered multiple values of dimensionality, sample size, correlation, regulation strength, and others. The adopted alternatives are possibly the most popular in data analysis and thus warrant direct comparison. In the analysis of TCGA data, the ANCut analysis results are more desirable. It is noted that the SA algorithm encourages but not forces clusters with comparable sizes. Examining the GO terms suggests that the clustering results can be biologically sensible.
This study inevitably has multiple limitations, including a lack of theoretical investigation, still limited simulations, and a lack of functional analysis of the data analysis results. More extensive analysis, numerical studies and comparisons will be conducted in the future. We note that the strategy of using regulator information to assist GE clustering is not limited to the NCut technique. In Additional file 1, we briefly discuss the possibility of trying such a strategy with Kmeans. We focus on the ANCut in this study and defer systematic development of the assisted clustering analysis to future research.
Conclusions
In this study, a new assisted analysis strategy is developed for clustering gene expression data. Advancing from the existing methods that focus on gene expression data only, information of regulators is utilized in clustering gene expressions. The proposed method carries a wealth of information and can reveal more accurate clustering. Experiment results on simulation and two TCGA datasets show the competitive performance of the proposed method with respect to the alternatives. This study provides a new venue for analyzing gene expression data. It is noted that although our analysis focuses on clustering of GEs, the proposed method is potentially applicable to other types of genetic measurements, such as proteins with their regulators.
Abbreviations
 ANCut:

Assisted normalized cut
 AugK:

Augmented Kmeans
 CNA:

Copy number alteration
 GE:

Gene expression
 GO:

Gene ontology
 NCut:

Normalized cut
 SA:

Simulated annealing
 Spec:

Spectral clustering
 SKCM:

Skin cutaneous melanoma
 TCGA:

The cancer genome atlas
References
Pers TH, Karjalainen JM, Chan Y, Westra HJ, Wood AR, Yang J, Lui JC, Vedantam S, Gustafsson S, Esko T, et al. Biological interpretation of genomewide association studies using predicted gene functions. Nat Commun. 2015; 6:5890.
Mar JC, Wells CA, Quackenbush J. Defining an informativeness metric for clustering gene expression data. Bioinformatics. 2011; 27(8):1094–100.
Napoleon D, Pavalakodi S. A new method for dimensionality reduction using kmeans clustering algorithm for high dimensional data set. Int J Comput Appl. 2011; 13(7):41–6.
Ding C, Li T. Adaptive dimension reduction using discriminant analysis and kmeans clustering. In: Proceedings of the 24th International Conference on Machine Learning. ACM: 2007. p. 521–8.
Kerr G, Ruskin HJ, Crane M, Doolan P. Techniques for clustering gene expression data. Comput Biol Med. 2008; 38(3):283–93.
Wang YR, Huang H. Review on statistical methods for gene network reconstruction using expression data. J Theor Biol. 2014; 362:53–61.
Wiwie C, Baumbach J, Röttger R. Comparing the performance of biomedical clustering methods. Nat Methods. 2015; 12(11):1033–8.
Ma S, Huang J. Regularized gene selection in cancer microarray metaanalysis. BMC Bioinforma. 2009; 10:1.
Hanisch D, Zien A, Zimmer R, Lengauer T. Coclustering of biological networks and gene expression data. Bioinformatics. 2002; 18(suppl 1):145–54.
Xu C, Liu Y, Wang P, Fan W, Rue TC, Upton MP, Houck JR, Lohavanichbutr P, Doody DR, Futran ND, et al. Integrative analysis of DNA copy number and gene expression in metastatic oral squamous cell carcinoma identifies genes associated with poor survival. Mol Cancer. 2010; 9(1):143.
Sun Z, Asmann YW, Kalari KR, Bot B, EckelPassow JE, Baker TR, Carr JM, Khrebtukova I, Luo S, Zhang L, et al. Integrated analysis of gene expression, CpG island methylation, and gene copy number in breast cancer cells by deep sequencing. PloS One. 2011; 6(2):17490.
Cutts RJ, Ullah AZD, Sangaralingam A, Gadaleta E, Lemoine NR, Chelala C. Ominer: an integrative platform for automated analysis and mining ofomics data. Nucleic Acids Res. 2012; 40(W1):560–8.
Zhu R, Zhao Q, Zhao H, Ma S. Integrating multidimensional omics data for cancer outcome. Biostatistics. 2016; 17(4):605–18.
Shi X, Zhao Q, Huang J, Xie Y, Ma S. Deciphering the associations between gene expression and copy number alteration using a sparse double laplacian shrinkage approach. Bioinformatics. 2015; 31(24):3977–83.
Zhou Y, Wang P, Wang X, Zhu J, Song PXK. Sparse multivariate factor analysis regression models and its applications to integrative genomics analysis. Genet Epidemiol. 2017; 41(1):70–80.
Zou H, Hastie T. Regularization and variable selection via elastic net. J R Stat Soc Ser B. 2005; 67:301–20.
Jiang Y, Shi X, Zhao Q, Krauthammer M, Rothberg B, Ma S. Integrated analysis of multidimensional omics data on cutaneous melanoma prognosis. Genomics. 2016; 107(6):223–30.
Swain PS, Elowitz MB, Siggia ED. Intrinsic and extrinsic contributions to stochasticity in gene expression. Proc Natl Acad Sci. 2002; 99(20):12795–800.
Shi J, Malik J. Normalized cuts and image segmentation. IEEE Trans Pattern Anal Machine Intell. 2000; 22(8):888–905.
Xing EP, Karp RM. Cliff: clustering of highdimensional microarray data via iterative feature filtering using normalized cuts. Bioinformatics. 2001; 17(suppl 1):306–15.
Yu Z, Luo P, You J, Wong HS, Leung H, Wu S, Zhang J, Han G. Incremental semisupervised clustering ensemble for high dimensional data clustering. EEE Trans Knowl Data Eng. 2016; 28(3):701–14.
Singh A, Soltani M. Quantifying intrinsic and extrinsic variability in stochastic gene expression models. Plos One. 2013; 8(12):84301.
Raser JM, O’shea EK. Noise in gene expression: origins, consequences, and control. Science. 2005; 309(5743):2010–13.
Bertsimas D, Tsitsiklis J, et al. Simulated annealing. Stat Sci. 1993; 8(1):10–5.
Meinshausen N, Bühlmann P. Stability selection. J R Stat Soc Ser B. 2010; 72(4):417–73.
Ng AY, Jordan MI, Weiss Y, et al. On spectral clustering: Analysis and an algorithm. Adv in Neural Information Processing Systems. 2002; 2:849–56.
Quackenbush J. Computational analysis of microarray data. Nat Rev Genet. 2001; 2(6):418–27.
Grün D, Lyubimova A, Kester L, Wiebrands K, Basak O, Sasaki N, Clevers H, van Oudenaarden A. Singlecell messenger RNA sequencing reveals rare intestinal cell types. Nature. 2015; 525(7568):251–5.
Stampfel G, Kazmar T, Frank O, Wienerroither S, Reiter F, Stark A. Transcriptional regulators form diverse groups with contextdependent regulatory functions. Nature. 2015; 528(7580):147–51.
McLendon R, Friedman A, Bigner D, Van Meir EG, Brat DJ, Mastrogianakis GM, Olson JJ, Mikkelsen T, Lehman N, Aldape K, et al. Comprehensive genomic characterization defines human glioblastoma genes and core pathways. Nature. 2008; 455(7216):1061–8.
Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, Davis AP, Dolinski K, Dwight SS, Eppig JT, et al. Gene ontology: tool for the unification of biology. Nat Genet. 2000; 25(1):25–9.
Consortium GO, et al. Gene ontology consortium: going forward. Nucleic Acids Res. 2015; 43(D1):1049–56.
Tibshirani R, Walther G, Hastie T. Estimating the number of clusters in a data set via the gap statistic. J R Stat Soc Ser B. 2001; 63(2):411–23.
Acknowledgements
Not applicable.
Funding
This study has been partly supported by NIH awards R21CA191383, R01CA204120, and P50CA121974, and NSFC award 61402276. National Bureau of Statistics of China award 2016LD01.
Availability of data and materials
The TCGA datasets used in this study are available at website https://tcgadata.nci.nih.gov/docs/publications/tcga/?.
Author information
Authors and Affiliations
Contributions
All authors were involved in methodological development. SH conducted programming and numerical study. All authors read and approved the final manuscript.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Additional file
Additional file 1
This file contains additional discussions and numerical results. (PDF 122 kb)
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Teran Hidalgo, S.J., Wu, M. & Ma, S. Assisted clustering of gene expression data using ANCut. BMC Genomics 18, 623 (2017). https://doi.org/10.1186/s1286401739901
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s1286401739901