BACKGROUND
systems biology is a descriptive paradigm in which one of the main concerns is how genes work together to form subsystems.
a basic assumption within this context is that genes which are co-expressed are often in the same subsystem  <cit> .
gene co-expression networks  are graph-based models used to express such subsystems.
construction of these networks is usually based on co-variation in expression within groups of genes across samples  <cit> .
they are graphs in which nodes are genes and edges represent interactions between them.
typically, the edges are undirected, in the sense that causality  is unassigned.
edges may be weighted and/or signed, thus indicating the strength of relationship between pairs of genes and up/down regulated interactions depending on the sign.
topological considerations, such as the number or relevance of connections for each node, can distinguish some nodes as highly interconnected  and central nodes within the system being modelled.
gcn can be used to make in silico functional predictions about genes.
the guilt by association  paradigm  <cit>  is used to predict function for genes that are not sufficiently studied and annotated using gcns.
gba assumes that genes that strongly co-express must share functionality, thus we can use well-characterised genes to assign function to those that are not.
groups of genes that tightly co-express are usually seen as a single functional unit.
on this basis, in the same way that single genes are used in association mapping with phenotype, convenient mathematical representations of groups of genes can be useful for multi-gene association mapping with phenotype as well.
one of the most widely used pipelines for gcn construction is weighted gene co-expression network analysis  .
it works in two main steps.
in the first step it constructs a network n of gene-gene co-expression in the form of a squared n×n matrix, where n is the number of genes in the study and each n is the interaction strength between the corresponding pair of genes .
in the second step, this matrix is used as the basis for obtaining a new squared distance matrix with the distance between genes, ready to be used for obtaining clusters.
and then such clusters can be used for multi-gene association mapping with traits or different downstream analyses .
this pipeline has been widely used and generated many fruitful insights into how genes interact within specific conditions .
in this paper we propose an improvement to the standard wgcna pipeline by a refinement of how the clusters  are generated.
this refinement is enabled through a hybrid clustering algorithm.
it uses the output of the conventional wgcna clustering as subsequent input to a k-means  <cit>  clustering algorithm for further refinement.
we will show that this hybrid scheme improves many interesting module properties paving the way to more accurate and potentially useful wgcna co-expression network analyses.
wgcna’s standard configuration uses hierarchical clustering .
in hc, a strong point is that the dendrogram structure eases the problem of finding a good number of clusters, k. moreover, the developers of wgcna include in the software an automated method to generate the appropriate number of clusters  <cit> .
on the other hand, a weak point of hc is that final results strongly depend upon how distances between clusters are compared.
furthermore, once the decision on which branch of the dendrogram a gene belongs to, this cannot be undone.
regarding k-means, a weakness in it is that the value of k  must be set prior to running the algorithm.
although there are techniques for setting it automatically, most of these are based on multiple random initialisations of centroids , so k is usually set arbitrarily.
it needs an initialisation of the centroids to start running.
a centroid is defined as an average representative of all the genes/points within the cluster such that all genes/points belonging to the cluster show minimum distance to that centroid in comparison to the other modules.
how we initialize these centroid will have a critical effect on the final result.
on the upside, k-means will search for the best centroids quickly and will quickly converge to an equilibrium situation .
the hybrid scheme we propose exploits the upsides from both approaches while alleviating their respective drawbacks.
k-means will move genes between modules thus effectively undoing premature decisions made by hc when assigning genes to sub-dendrograms.
we set the value of k equal to the number of modules discovered by hc and we initialise the centroids to the eigengenes generated by wgcna, thus taking advantage of hc to carry out sensible initialization .
implementation
the standard wgcna procedure
consider a gene expression profile matrix g
n×m where n is the number of samples for a given condition, m is the number of transcripts and each g in g gives the quantification of the j-th transcript within the i-th sample.
the standard wgcna  <cit>  procedure generates a squared adjacency matrix, between genes, based on their correlation.
depending on whether the adjacency is signed  or unsigned  we will obtain networks either reflecting the direction of co-regulation  or ignoring it, respectively.
adjacency is defined as a
d
j=|c
o
r|β for genes i and j. the β parameter is an integer that modulates how smooth is the transition between the lowest to the highest possible co-regulation between genes.
the wgcna methodology enables choosing β in such a way that the network shows a scale free topology  property  <cit>  .
this feature is commonly observed in biological networks.
from the adjacency values, a new matrix with the same dimensions is created, the topological overlap matrix .
this step alleviates the effect of noisy genes when obtaining the adjacency from correlation.
once the network is built through the tom, it is converted to a distance matrix  to use it as the basis for clustering .
a dynamic tree-cutting algorithm  <cit>  is then applied to the dendrogram to generate a partition p={p
<dig> …,p
k} of disjunct sets of genes.
thus, wgcna generates two main components which are useful for subsequent downstream analyses.
on the one hand, the tom gives, for the j-th row/column, the level of co-expression of gene j with all of the genes in the network.
the higher the value for a given  pair, the tighter the interaction between them.
furthermore, the sum of all row or column values for a gene, will give a measure of its overall level of co-regulation within the experimental condition, i.e.
its ‘hubbiness’.
thus, the tom is, in effect, the gcn.
the other component produced by wgcna is the partition of gene sets, p, created from the tom.
these partitions or modules often reflect cell types, common cellular functions or other biological subsystems reflecting, for example, immune function, or function related to the tissue under study  <cit> .
but the main utility of modules is to allow mapping gene groups to traits, when available.
following the wgcna standard methodology, this is performed by looking for significant correlations between traits and the module ‘eigengene’.
the eigengene summarizes the overall module activity in a given sample, and is obtained as the 1st pc component of the gene expression of genes belonging to the module.
improvement of hierarchical clustering with k-means
hc provides a convenient graphical representation of groupings that can be validated by biologists.
one can readily obtain a suitable number of clusters from such an approach by ‘cutting’ the dendrogram at different heights, either manually or via various automatic algorithms  <cit> .
but, as we explained above, the final dendrogram strongly depends on how we measure distance between clusters .
furthermore, once a gene falls under a subdendrogram, this decision cannot be modified under hc.
if we consider how wgcna manages modules and eigengenes, it is assumed that each gene is highly correlated with other genes in its module.
in other words, the module membership  of the gene in its own module, measured as the pearson correlation between its expression and the module eigengene, should be higher than it is for any other module.
however, we show here, from our real-data analyses, that 25% of genes would be better off in other modules .
in this paper, we propose a post-processing step based on k-means to overcome all these limitations.
it works on the partition p, leaving the tom unmodified.
the k-means algorithm  <cit>  is well known and works on the n dimensional sample space of m points  in an iterative fashion.
it starts by setting a value for k, the number of clusters to discover and k centroids, one for each cluster.
centroids are the representatives of each cluster, in such a way that a point  g belongs to cluster i if the distance of such point to the cluster centroid is the minimum among all distances to all k cluster centroids.
in standard k-means, given a partition of k modules, the the centroid for the i-th module c
i={c
i, <dig> …,c
i,n} is generated as follows
<dig> ci=1n∑j=1mgj,wheregi∈pi.
however, in wgcna, the notion of a centroid is substituted by that of an eigengene.
accordingly, our definition of k-means will use eigengenes as centroids.
the concept of distance is a central element of k-means.
it is important to note that distance in k-means is always defined between a point in the dataset  and a centroid .
euclidean distance is the the most commonly used distance in conventional k-means.
however, given that we are constructing co-expression networks based on correlations, distance cannot be euclidean.
modules should represent co-expressed genes  instead.
thus, and depending on the wgcna type of network, we should apply a distance between gene and eigengene based on the co-expression measure used.
we will limit our discussion to signed networks.
these specific types of networks will separate up- from down-regulated genes in different modules, which is usually of biological interest.
they are also convenient for downstream analysis as correlation of genes and eigengenes will be positive, which eases a posteriori analyses.
in signed networks, wgcna uses
<dig> co=12),
as a normalised measure of co-expression between the expression profile of a gene g
i and a eigengene e
g
j, where by default c
o
r() is the pearson correlation coefficient.
accordingly, we use 1−c
o as distance.
it is worth noting that hc needs a distance matrix between all genes, i.e.
1−t
o
m. k-means needs instead a computable distance definition between gene and eigengene.
finally, on the basis of this definition of centroid and distance, genes are reassigned to the partitions induced by the new centroids, iteratively.
if a stopping criterion is met, the algorithm finishes.
otherwise, a new iteration is performed.
we note that wgcna is computationally optimized to use pearson correlation.
other correlation measures are in principle possible, including spearman’s rank correlation coefficient.
however, in our own investigations we observe an increase in computation time of at least × <dig>  when using spearman correlation without seeing any conclusive improvement with respect to the biology of networks .
thus, throughout this paper we perform analyses using pearson correlation.
we propose a general procedure which obtains, from a g
n×m matrix of gene expression profiles from n samples and m genes, a clustering partition p of such genes by incorporating the standard wgcna process together with a post-processing of the partition obtained from it.
the original contribution of this paper is described in steps from  <dig> to  <dig> below.
step 1: initialization.
let g
n×m be a dataset of n samples and m genes for a given condition.
let d() denote a distance function between a gene in g and an eigengene.
let f
c() denote a function which takes a clustering partition p={p
<dig> …,p
k} as an argument and generates centroids
step 2: β=wgcna::picksoftthreshold$powerestimate
step 3: obtain a tom, given g and β
step 4: generate a partition p
hc={p
<dig> …,p
k} with 1−t
o
m as a distance matrix and with average linkage hierarchical clustering and dynamic cutting height.
step 5: let c={c
<dig> …,c
k} be a set of k vectors of n components which denote the centroids of the k-means clustering.
step 6: initialize c with f
c
step 7: create a new partition p
km={p
<dig> …,p
k} by assigning each gene g
i, 1≤i≤m to a p
j∈p
km such that d≤d, 1≤t≤k holds.
step 8: if the termination criterion holds, then stop.
if not, generate new centroids c with f
c and go to step  <dig>
note that, in this algorithm, we left d() and f
c() undefined.
however, within this paper we define d() according to eq.
<dig> and f
c() as the module eigengene.
computational complexity of the proposed approach
conventional wgcna gcn construction needs three sequential steps:  obtaining the soft threshold  to account for scale free topology, which has a computational complexity that depends on the number of genes and samples,  obtaining the tom matrix, which has a complexity o, where n is the number of genes, as it has to construct a n×n squared matrix of adjacencies between all n genes, and  hierarchical clustering, which has a complexity o. overall, wgcna’s computational complexity is o. wgnca’s space complexity is also o because it needs to maintain the tom in memory for hc to get the clusters.
the computational complexity of k-means fits well with wgcna’s complexity.
its time complexity is o where n is the number of genes, k is the number of clusters and it is the number of iterations.
assuming that k,i
t< <dig>  using it as a post-process is very affordable in terms of computation time.
note that k-means does not require the tom matrix in memory as the only distances it requires are between genes and eigengenes, and these we obtain on the fly by using eq.
<dig>
stopping criterion
it is reasonable to assume that a sufficiently high number of k-means iterations will always be able to decrease the number of misplaced genes  to  <dig>  on the other hand, the algorithm’s time complexity  means that it is possible to run a single k-means in a conventional laptop in a matter of a few minutes.
this means that we could, in principle, design a stopping criterion based on the minimun number of misplaced genes being set to  <dig>  however, we note that a situation may exist where the algorithm may fall into an infinite loop without reaching the desired state .
thus, the stopping criterion we include in the software package km2gcn tries to reach the desired value for misplaced genes but always within a limited number of iterations.
we did not observe the mentioned infinite loop situation in any of our experiments.
RESULTS
we wished to assess the ability of our method to define gene groups that genuinely reflect biological function.
this is non-trivial for the following reasons.
firstly, many genes are known to be pleiotropic, i.e.
a single gene can affect many traits  <cit> .
transcription factors are a good example of this  <cit>  but there are many other examples  <cit> .
by creating non-overlapping partitions we deliberately ignore this fact and implicitly assume a model in which genes are highly specialized .
secondly, we are limited by technology and sample availability from producing optimum estimates of gene expression profiles.
we therefore lack of all the necessary information to build the best model.
finally, if we wanted to evaluate the functional similarity of genes within a module, again we do not know all functions that all genes may play in any condition.
notwithstanding these caveats, we explored various approaches to provide a comprehensive and varied assessment of the effectiveness of our k-means hybrid method.
in “materials and methods for the gcns used for our evaluations” section we describe the datasets used in our investigations and the particular pipelines used to obtain the corresponding gncs.
in “dynamics of k-means when working on 1-tom distance space” section we show our hybrid approach  works.
in “is k-means doing a proper job?” section we digress to note that k-means actually optimizes the sum of squares of within cluster distance.
in “k-means improves the ‘eigengene’ as a tool for analysis” section we show that the proposed approach improves modules as a tool for mapping with traits.
in “k-means improves module preservation” section we suggest that k-means improves cluster similarity between conditions .
in “k-means detects more accurate partitions than wgcna in simulated data” section we compare the accuracy of k-means against wgcna on simulated data.
in “k-means improves functional enrichments” section we show that k-means improves a module’s functional characterization through well-known databases such as the gene ontology.
finally, in “k-means improves brain specific cell type marker enrichment” section we present results that suggest that gene markers for specific cell types show a better arrangement in partitions generated from k-means.
materials and methods for the gcns used for our evaluations
we evaluated gcns in two well-known datasets.
the first  is focused on brain tissue exclusively and it is based on affymetrix human exon v <dig> microarray expression profiles from  <dig> brain tissues.
this dataset is well suited for evaluating the k-means extension to wgcna because it is well known, it comprises  <dig> different brain regions and gcn networks created with the standard wgcna method have been published  <cit> .
the procedure used to create the gcns is as follows.
sample outliers were identified by visual inspection after clustering the samples using hierarchical clustering with euclidean distance as the distance measure.
the majority of the identified outliers had low interarray correlation, which is defined as the pearson correlation coefficient of the expression levels for a given pair of transcripts using all available data available .
after outlier removal, the same process was repeated to check for additional outliers.
the gcn constructed was of signed type, with β= <dig> for all tissues.
using these settings, the hc wgcna partition was created using  <dig>  transcripts  passing quality control.
once the partition was created,  <dig> additional transcripts  were assigned to modules based on their highest module membership.
each partition was refined afterwards with k-means.
the second dataset is gtex  <cit> , which is one of most comprehensive human datasets currently available for multi-tissue transcriptomics.
the gtex v <dig> gene expression dataset comprises  <dig>  samples unevenly distributed across  <dig> post-morten human tissues.
we created networks for  <dig> tissues.
in sequential steps, starting from rpkm  <cit>  values of gene-level quantification provided by gtex, we selected all tissues with more than  <dig> samples.
for each tissue we retained ensembl genes with r
p
k
m> <dig>  seen in more than 80% of the samples.
this produced a variable set of genes for each tissue, with a minimum of  <dig>  for skeletal muscle tissue and a maximum of  <dig>  for testis.
we applied batch, gender, age and rin as known covariates for data correction and to account for unknown covariates we applied sva  <cit>  axes.
for each dataset of filtered rpkm values, we applied the sva r package using svaseq() to generate sva axes.
for network construction we used the residuals obtained by regressing the rpkm expression values with the known and unknown covariates with a generalized linear model.
to construct the networks, we applied the algorithm introduced in “improvement of hierarchical clustering with k-means” section.
note that the differences between the ukbec and gtex networks are important and makes them well suited and complementary for the purpose of our study.
the ukbec gene expression dataset is microarray based, while the gtex gene expression dataset is based on rna-seq technology, with rpkm quantification.
ukbec networks are restricted to  <dig> brain tissues while gtex networks cover  <dig> tissues, including  <dig> brain tissues .
in summary, we have  <dig> gtex gncs,  <dig> brain specific ukbec gcns, gtex sample sizes range from n= <dig> to n= <dig> ; ukbec sample sizes range from n= <dig> to n= <dig> ; we have a variable number of genes used in the gtex gcns, in the range  <dig>  to  <dig>  for skeletal muscle and testis respectively ; and the same  <dig>  probes for all  <dig> ukbec gcns.
finally, note there is a much higher variability in the number of modules per gcn in gtex,  <cit>   than in ukbec,  <cit>  .
please note that throught the paper we use abbreviations to refer to tissues.
please see table  <dig> for the correspondence between abbreviations and brain region names.
dynamics of k-means when working on 1-tom distance space
as outlined in “improvement of hierarchical clustering with k-means” section, our proposed algorithm does not modify the distance matrix  but acts later, on the partition p={p
<dig> …,p
k} taking k from the number of modules discovered by the hc used within wgcna.
k-means acts iteratively creating centroids from modules and deciding for each gene, on the basis of the new centroids, which one is nearest to the gene.
if, in the current iteration, the gene is nearest to a different centroid, then the k-means algorithm assigns it to the corresponding module.
thus, in each iteration a new partition is generated with the changes applied to the former partition.
i∈p, on average 30% of other modules within the gcn contribute with genes to its final gene set configuration.
genes that leave their hc module have a module membership at that module of  <dig>  on average, with standard deviation of  <dig> .
genes arriving to a module for the final k-means partition show an average mm on arrival,  <dig> , with standard deviation  <dig> .
fig.
<dig>
upper plot shows the evolution of the number of moved genes  between any pair of modules p
i and p
j across k-means iterations  for ukbec-microarray dataset.
bottom plot shows the average module membership of genes  moved  across iterations  for the ukbec-microarray dataset in comparison with average module membership for all the genes
the lower panel of fig.
<dig> focuses on the  <dig> ukbec gcns, and on how mm evolves with iterations.
dashed lines show, for each tissue, the average mm of moved genes, at each iteration, defining mm with reference to the original wgcna partition.
initially, the algorithm focuses on moving genes with very low mm, but following this it focuses on genes with higher mm and then stabilizes.
the solid lines show the average mm of all the genes in the network for each iteration.
this dramatically increases over the first iterations and then smoothly and monotonically increases across additional iterations.
this suggests that, over time, genes’ assignment to a module becomes stronger.
is k-means doing a proper job?
k-means is designed to optimize the sum of squares of within clusters distance  <cit> .
we define the within cluster distance, denoted with w, for a partition p as
<dig> w=∑k=1k∑p=k||gi−ck|| <dig>
where k is the number of clusters within p, p refers to the cluster that g
i belongs to, and c
k is the eigengene for the k-th cluster.
we could, alternatively, define the distance between clusters, b. note that for a given set of genes, we can obtain the sum of distances between all gene pairs.
if we denote this measure by t for a given gene expression profile g, we can decompose it into t=b+w for any given p obtained from g .
this means that maximizing the between-clusters distance is equivalent to minimizing the within-cluster distance.
we observe that k-means monotonically decreases w  increases) across iterations .
the algorithm generates a higher w at the early iterations, which decreases to a lower level in later iterations.
this behaviour is in line with the shape of the gene changes curves of the upper plot in fig.
<dig>  higher number of moved genes imply higher decreasing rate of within cluster distance.
fig.
<dig> the within cluster distance evolution during the k-means runs for the ukbec datasets
this behaviour is also in accordance with what we see in fig.
<dig>  this plot shows, for each module p
i∈p, and the specific case of ukbec’s cerebellum cgn, the distance between the eigengenes for the same module, as they are created during successive iterations.
over time, the eigengene vectors stabilize across iterations suggesting that cluster definition becomes stable.
fig.
<dig> euclidean distance of successive module eigengenes along the k-means iterations for cerebellum samples for ukbec datasets
k-means improves the ‘eigengene’ as a tool for analysis
one of the main applications of wgcna partitions is searching for associations between gene clusters and traits.
traits are usually given as a vector of n components where n is the number of samples.
on the other hand, modules are comprised of m vectors of the same form where each vector is the expression profile for the corresponding gene.
to assess the correlation between the trait and the gene module, wgcna transforms module p
i into an eigengene .
from this, the correlation between the trait and the eigengene can be easily obtained.
there are several applications for the eigengene.
for example, it can be used to provide a measure of how strong is the membership of each gene g∈p
i to the i-th module, by correlating its expression with the eigengene, resulting in the mm of g. let this module membership be denoted with m for gene g and module i. it is assumed that a good p would be one such that, the number of genes g with m<m when g∈p
i, for any i≠j, is low.
let us call such genes ‘misplaced genes’.
we would prefer a partition in which the number of misplaced genes is minimum.
to assess the number of misplaced genes after applying k-means, we performed investigations in the  <dig> ukbec and  <dig> gtex tissues.
in ukbec, using k-means with only  <dig> iterations , we get a maximum of  <dig> misplaced genes in the putamen and a minimum of  <dig> for occipital cortex and a mean of  <dig> misplaced genes per partition.
in the wgcna partitions, the maximum number of misplaced genes is  <dig> in temporal cortex, the minimum is  <dig> in white matter and an average of  <dig> genes;  <dig> times more than with the k-means algorithm.
in gtex clustering partition modules, the average number of misplaced genes in modules from a wgcna partition is  <dig>  after applying k-means, it is only  <dig> .
k-means improves module preservation
one component of wgcna provides a convenient tool for the analysis of module preservation  <cit> .
given a partition, p, constructed from a network obtained from a given set of samples s, we can test whether the features of each module p
i∈p  are preserved in an alternative set of samples s
′ .
preservation analysis is based on estimating, for some statistic of interest, differences between what is observed and what is obtained by random permutation.
for example, one statistic of interest is the gene correlation with the eigengene .
through a simple transformation one can check whether the values obtained in the reference network are maintained  for the same genes within the other network.
wgcna uses the ‘z-summary’ statistic as a general summary of all the different statistics used.
to assess the effect of k-means on z-summary, we performed the same investigations on both  <dig> ukbec brain tissues and on the  <dig> gtex brain tissues.
note we focus on brain tissues within gtex as comparison of preservation only makes sense for tissues that are similar.
within each ukbec and gtex tissue gcn, we compared the preservation of all the partitions generated by wgcna alone with the preservation obtained by applying k-means to each of them.
a permutation analysis on  <dig> tissues generates, for each tissue t, and for each module p
i within the corresponding tissue network, a vector of  <dig> z-summary values corresponding to the preservation of p
i in the other  <dig> tissues.
table  <dig> number of new modules from a tissue  that are preserved on another tissue  after applying the k-means to the standard wgcna partitions
from table  <dig> it is apparent that there is an overall increase in the number of modules preserved under k-means.
in the ukbec gncs, there is an improvement in  <dig> cases , no improvement in  <dig> cases , and only case with a worse preservation .
the average improvement in modules preserved for ukbec is  <dig> .
in the gtex gcns, there is an improvement in module preservation in  <dig> cases , no improvement in only  <dig> cases  and a decreased preservation in just  <dig> cases.
the average number of modules improved by the k-means method is  <dig>  .
this suggests that k-means creates less noisy modules as similarities between tissues are more apparent.
finally, it is worth noting that each tissue is expected to have specific modules, i.e.
modules that will be poorly preserved in other tissues because they are exclusive from that tissue, reflecting study-specific or sample-specific gene subsystems.
k-means detects more accurate partitions than wgcna in simulated data
we wanted to test whether k-means improves the accuracy of partitions p with respect to those obtained under standard wgcna.
to this end we investigated networks based on ‘synthetic data’.
the wgcna package provides a gene expression profiling simulation method simulatedatexp(), which is a convenient method for generating artificial data sets that mimic the properties of real datasets.
the simulation method works with the eigengene of gene expression for each gene belonging to the module.
the simulation method requires, as arguments:  a matrix with the eigengene for each module;  the proportion of the total gene pool that one will find within each module; and  the number of genes to be simulated.
note that the number of samples we want to simulate appears implicitly as the length of each eigengene .
the method returns two elements:  a gene expression profiling data set, let us denote it with d, that we can use to construct gcns; and  the ideal clustering partition of the simulated gene expression profiling, here denoted by p. thus, if we rely on the effectiveness of this simulation method, then a simulated data set d, we will prefer a gcn construction algorithm a to algorithm b if the distance between a and p is smaller than between b and p, where a and b are the clustering partitions we get after constructing gcns on d with a and b, respectively.
the accuracy of an algorithm a is defined by the similarity of the theoretical optimal partition within the synthetic data to the partition constructed by a.
in order to test whether k-means performs any better than standard wgcna on simulated data, we constructed a plausible set of simulated gene expression profiles.
we used gtex and test with them both k-means and standard wgcna on gcn construction.
the accuracy of an algorithm a will be defined as how similar are the theoretical optimal partition within the synthetic data, and the partition constructed by a.
in order to test whether k-means performs any better than standard wgcna on simulated data, we constructed a plausible set of simulated gene expression profiles.
we used gtex standard wgcna gcns  as the simulation seed for the generation of a synthetic gene expression profile.
we focused on the gtex dataset rather than ukbec, because the  <dig> gcns comprise a usefully heterogeneous network dataset.
the simulated data process produced a gene expression profile and a theoretical ideal clustering partition for such profile.
we used this theoretical ideal partition to evaluate standard wgcna and k-means accuracy.
to estimate accuracy we use three different statistics:  the rand  <cit>  index, also implemented within wgcna, the jaccard coefficient and the similarity index  <cit> , all of them implemented within clv r package.
results for all the experiments appear in additional file  <dig>  each row corresponds to a gtex tissue, the randsimvswgcna column corresponds to the rand index between the ideal partition and that obtained with wgcna on the simulated data.
the randsimvskm column corresponds to the same index when using k-means.
the other four columns correspond to the jaccard coefficient and the similarity index.
the k-means refinement generate higher values in all the cases for all three indexes.
these results are illustrated in fig.
<dig>
fig.
<dig> results on performance of standard wgcna and k-means on  <dig> simulated data sets that used the gtex wgcna gncs as seed for simulation.
we display the same results using three different indexes of similarity between cluster partitions.
the k-means method outperforms standard wgcna with all three indexes used
k-means improves functional enrichments
the gene ontology  <cit>  is a curated database for gene annotation which can be used for the functional characterization of gene sets.
given a set of genes , and a subset of those genes , an enrichment analysis can be performed on go annotations  <cit>  to search for terms in the ontology that are significantly enriched in the subset of genes relative to the full set.
the number and strength of significant terms obtained in this way can be used to measure the biological functionality of the module.
given two different partitions p and p
′ created from the same tom, we would prefer the partition that generates more significant go terms if we assume go to reflect a biological ground truth as this would suggest that the preferred partition makes more biological sense.
we used the gprofiler r package  <cit>  to obtain enrichment p-values, avoiding eia  terms in go and requiring a correction for multiple testing with gscs, as developed by the authors of the package.
we describe below a series of investigations to characterize the improvement in a module’s biological functionality.
global annotation term significance
consider a partition p={p
<dig> …,p
k} of genes arranged into modules p
i, 1≤i≤k.
now suppose we want to perform a gene set enrichment analysis on each p
i∈p based on the gene ontology.
go is a list of ontological terms, organised into three main branches: bp , mf  and cc .
genes within the database will be associated with a number of terms from each branch.
thus, for each term in go, and given the list of genes in p
i, we can apply a contingency test, e.g.
fisher exact test  <cit> , under the null hypothesis that the genes in p
i show no significant overlap with the set of genes associated with the term.
with an appropriate correction for multiple testing, we define as significant the association of the list of genes in p
i with the corresponding term, when the corrected p-value is < <dig> .
we then aggregate all these p-values for a module in a single measure of significance as follows.
for each p
i∈p, we use
<dig> sgo=∑pvaluej∈test−log <dig>
where t
e
s
t is the set of p-values, p
v
a
l
u
e
j, of significant terms associated with the genes in partition p
i, emerging from the analysis.
in this way,
<dig> sgo=∑pi∈psgo
can be used to aggregate all the biological signals  of a whole partition p. given a choice of partitions, we prefer p to p
′ when s
go>s
go.
a displays for each ukbec and gtex gcn, the relative improvement between the standard wgcna partition, p and the k-means partition, p
′ by
sgosgo− <dig>
fig.
<dig> the left plot’s light blue blue bars show the percentage of relative improvement by k-means with respect to wgcna s
go statistic.
values in red  are those that k-means fails to improve.
the right plot shows cell type enrichment improvement in the same way, for the  <dig> ukbec gcns and the  <dig> gtex brain networks.
again, values in red are those that k-means fails to improve
the average improvement is 13% .
overall, there is improvement in all ukbec tissues and in  <dig> out of  <dig> gtex gcns, and the overall improvement is significant .
does a higher enrichment implies less informative modules?
high values of the s
go index are of interest, as we prefer a partition p over p
′ if s
go>s
go.
however, it is possible that modules show better s
go values after k-means because the module have more annotation terms that are generic, and therefore less descriptive about the specifics of the tissue studied.
in order to assess this, we applied the notion of information content  <cit> .
we used the gosim package  <cit>  which applies information-based metrics to gene ontology terms.
the metric i
c for a term t belonging to an ontology is defined as:
<dig> ic=−logp,
where p is the probability of observing t within the annotations available within that ontology.
ideally, we prefer modules with more go terms, which are more significant  and more informative .
from previous sections we know we have more significant networks thanks to s
go.
but is k-means capable of not only improving significance but also of maintaining the level of information of the modules if not increasing it?
fig.
<dig> relation between frequency of appearance of go annotation terms across all gtex gcns and ic .
terms appearing more times tend to have lower ic.
regression lines show that k-means gets better ic values for highly repetitive terms
is the overall ic obtained by k-means degraded as a consequence of obtaining more significant terms per gcn in comparison with standard wgcna?
to assess this we regressed the information content of the significant annotation terms against the frequency of appearance in the gcn annotation sets.
we found a tendency towards higher ic in k-means gcns.
this suggests that k-means annotations are more specific, and therefore more useful.
is the increase in enrichment better than random?
in “dynamics of k-means when working on 1-tom distance space” section we noted that one of the changes within cluster partitions after applying k-means is that module sizes change and many modules will increase their size considerably.
it is fair to assume that modules increasing their size in genes, will also increase their s
go enrichment.
there is a significant pearson correlation between increase in module size after k-means and increase in number of significant annotation terms .
the question arises, therefore, of what is the real contribution of k-means in comparison to a random shift of genes between modules?
in order to answer this question, for each of the  <dig> gtex tissues and their corresponding wgcna and k-means partitions, we identified those genes that were changed at the wgcna partition to create the k-means one.
then, in a single step, we randomly assigned these genes to other modules in such a way that we kept the same module sizes obtained with k-means.
via this algorithm, we produced new partitions in which the genes that remained unchanged from wgcna to k-means stayed at the same modules, but those genes that were changed by k-means were again changed but this time in a random fashion.
go statistic.
plot  shows the number of significant go terms.
fig.
<dig> effect of random assignment of genes selected by k-means, on a wgcna partition, to be changed from one module to another.
plot  refers to s
go values and  to number of significant terms
in 89% of the modules, k-means finds the same number  or more  significant go terms than the random placement of misplaced genes .
88% of the final modules show equal  or better  s
go index using k-means compared to random .
aggregating the results by tissue, k-means placement performs better in all the cases.
interestingly, the random placement of misplaced genes prevents enrichment at the wgcna partition it starts with .
this is important because even though many genes at the wgcna partition are not touched by the random approach, moving genes randomly will nevertheless worsen these genes’ functional annotations.
this suggests that both the number of significant terms and the s
go index have a reasonable sensitivity.
k-means improves brain specific cell type marker enrichment
one interesting property of wgcna gcns is that partitions created from them can be useful when studying cell-type specific gene networks.
in studies where samples come from bulk tissue, it is most likely that these samples will be comprised of different cell material.
in consequence, the gene expression profiles obtained from them should reflect this heterogeneity in some way.
wgcna’s gcns handle this heterogeneity in an elegant and convenient manner: they often generate gene clusters within partitions which are specialized on a given cell type, i.e.
they present a highly significant enrichment of markers  for a given cell type  <cit> .
we wanted to assess the effects of k-means on this particular feature.
to do this, we used three different resources defining cell-type specific gene sets.
these were wgcna’s brain lists,  <cit> , and two alternative brain specific sources, labelled here external  <cit>  and cahoy  <cit> .
we evaluated each partition’s modules from the  <dig> ukbec gcns and the  <dig> brain tissue gtex gcns, using both standard wgcna and k-means.
this evaluation generated two matrices of p-values , with each gene dataset in a row and each specific module from each of the  <dig> networks in a column.
p-values reflect a fisher’s exact test for whether there is significant concentration of the corresponding gene sets in the tested module.
we include in additional files  <dig> and  <dig>  the results for standard wgcna and k-means, on the  <dig> ukbec gcns.
note that, in these plots, both columns and rows have been clustered based on −l
o
g <dig> so it can be better seen how modules from different tissues cluster together at columns, and also how different gene sets cluster among rows.
these heat-maps reveal strongly clustered areas corresponding to groups of cell-type specific genes sets within most, if not all, of the tissues.
more specifically, we see four groups cell-type specific gene sets corresponding to microglia, astrocytes, oligodendrocytes and neurons .
in the ukbec k-means heat-map , using a significance cut-off of 10− <dig> , almost 65% of the modules show cell type enrichment .
within these,  <dig> modules show a single cell type signal.
in the wgcna heat-map , 63% of the modules show cell type enrichment , with  <dig> showing single cell type signals.
b compares the two enrichment matrices, by aggregating all the cell-type enrichment −l
o
g
<dig> transformed p-values as we did for the gene ontology enrichment in “k-means improves functional enrichments” section.
each bar represents the sum of all values of the corresponding heat-map, for modules of the given tissue.
according to this, we always see an improvement in ukbec networks and in most of the the gtex networks.
the overall improvement is significant .
CONCLUSIONS
our study shows that an additional k-means step, when used as an adjunct to wgcna, improves the partitions generated from gene co-expression networks.
our method is not an alternative to wgcna, instead it is an additional step to the standard wgcna pipeline.
indeed, our method can be applied to any general hierarchical clustering algorithm, and as such it could be usefully applied to any hierarchical clustering based approach for network generation, not just gene co-expression networks.
we evaluated our method using two contrasting gene expression datasets representing a variety of different tissues, one obtained with microarray technology , and the other with rna-seq .
using a variety of approaches, we demonstrate improved performance of our k-means method in both datasets.
furthermore, we also demonstrate improved performance using simulated data generated from the gtex dataset.
we show via these analyses that it is possible to obtain better partitions for the same networks via our k-means method.
our method generates modules with fewer misplaced genes with respect to their eigengene, and this implies that the eigengene is a better representative of the phenomena hidden behind the particular set of genes belonging to the module.
using gene ontology enrichment analyses, we also show that our partitions are enriched for biological functionality.
statistically significant s
go enrichment is seen in all  <dig> ukbec cgns and in  <dig> out of the  <dig> gtex gcns.
our partitions have improved modules preservation, which also suggests that the clustering is more accurate from a biological point of view.
although some gene modules are specific of each tissue , it is a desirable property of most gcn partitions to be highly replicable under the assumption that a preserved module is more likely to be a genuine module.
our analyses suggest that k-means favours the creation of more genuine modules and these results are seen in both ukbec and gtex gcns.
our k-means method also creates partitions in which gene sets representative of specific brain cell types are seen in modules with increased statistical significance.
this suggests, once again, more biologically genuine modules.
gcn construction is likely to become an increasingly important analysis, as genomics and transcriptomics are increasingly applied to aid clinical diagnosis and prognosis.
methods that generate more reliable and robust gene gcns will enable improved prediction of inter-gene relationships and gene function, with a variety of applications.
availability and requirements
ukbec data  <cit>  has accession code gse <dig>  all information about tissues, samples and quality control can be found there.
gtex rpkm gene expression v <dig> was used in this paper and downloaded from the gtex portal: http://gtexportal.org/home.
regarding the software we present here, this is the availability and requirements.
project name: km2gcn
project homepage:
https://github.com/juanbot/km2gcn
operating system: linux/windows/mac
programming language: r
other requirements: wgcna r package and gprofiler
license: lgpl
additional files
additional file  <dig> lists tissues, samples and genes used for the creation of each gcn.
additional file  <dig> includes comparative results of standard wgcna and k-means on simulated data.
additional file  <dig> the second one corresponds to k-means.
both on ukbec datasets.
values higher than  <dig> are set to  <dig>  colors at the top of columns correspond to tissues, the tissue legend is at the bottom of columns and cell marker gene set used on the right side.
additional file  <dig> heat-maps showing −l
o
g
<dig> from ficher’s exact test on significant concentration of specific cell marker gene sets  on each tissue module .
the one within the additional file  <dig> corresponds to the standard wgcna.
abbreviations
bpbiological process category of the gene ontology
cccellular component category of the gene ontology
eiaelectronic inferred annotations
gbaguilty by association
gcngene co-expression network
gogene ontology
gtexgenotype-tissue expression project
hchierarchical clustering
mfmolecular function category of the gene ontology
mmmodule membership of a gene
sftscale free topology
tomtopological overlap matrix
ukbecunited kingdom brain expression consortium
wgcnaweighted gene-co expression network

 systems biology is a descriptive paradigm in which one of the main concerns is how genes work together to form subsystems. a basic assumption within this context is that genes which are co-expressed are often in the same subsystem. gene co-expression networks are graph-based models used to express such subsystems. construction of these networks is usually based on co-variation in expression within groups of genes across samples. they are graphs in which nodes are genes and edges represent interactions between them. typically, the edges are undirected, in the sense that causality is unassigned. edges may be weighted and/or signed, thus indicating the strength of relationship between pairs of genes and up/down regulated interactions depending on the sign. topological considerations, such as the number or relevance of connections for each node, can distinguish some nodes as highly interconnected and central nodes within the system being modelled. gcn can be used to make in silico functional predictions about genes. the guilt by association paradigm is used to predict function for genes that are not sufficiently studied and annotated using gcns. gba assumes that genes that strongly co-express must share functionality, thus we can use well-characterised genes to assign function to those that are not. groups of genes that tightly co-express are usually seen as a single functional unit. on this basis, in the same way that single genes are used in association mapping with phenotype, convenient mathematical representations of groups of genes can be useful for multi-gene association mapping with phenotype as well. one of the most widely used pipelines for gcn construction is weighted gene co-expression network analysis. it works in two main steps. in the first step it constructs a network n of gene-gene co-expression in the form of a squared n×n matrix, where n is the number of genes in the study and each n is the interaction strength between the corresponding pair of genes. in the second step, this matrix is used as the basis for obtaining a new squared distance matrix with the distance between genes, ready to be used for obtaining clusters. and then such clusters can be used for multi-gene association mapping with traits or different downstream analyses. this pipeline has been widely used and generated many fruitful insights into how genes interact within specific conditions. in this paper we propose an improvement to the standard wgcna pipeline by a refinement of how the clusters are generated. this refinement is enabled through a hybrid clustering algorithm. it uses the output of the conventional wgcna clustering as subsequent input to a k-means clustering algorithm for further refinement. we will show that this hybrid scheme improves many interesting module properties paving the way to more accurate and potentially useful wgcna co-expression network analyses. wgcna’s standard configuration uses hierarchical clustering. in hc, a strong point is that the dendrogram structure eases the problem of finding a good number of clusters, k. moreover, the developers of wgcna include in the software an automated method to generate the appropriate number of clusters. on the other hand, a weak point of hc is that final results strongly depend upon how distances between clusters are compared. furthermore, once the decision on which branch of the dendrogram a gene belongs to, this cannot be undone. regarding k-means, a weakness in it is that the value of k must be set prior to running the algorithm. although there are techniques for setting it automatically, most of these are based on multiple random initialisations of centroids, so k is usually set arbitrarily. it needs an initialisation of the centroids to start running. a centroid is defined as an average representative of all the genes/points within the cluster such that all genes/points belonging to the cluster show minimum distance to that centroid in comparison to the other modules. how we initialize these centroid will have a critical effect on the final result. on the upside, k-means will search for the best centroids quickly and will quickly converge to an equilibrium situation. the hybrid scheme we propose exploits the upsides from both approaches while alleviating their respective drawbacks. k-means will move genes between modules thus effectively undoing premature decisions made by hc when assigning genes to sub-dendrograms. we set the value of k equal to the number of modules discovered by hc and we initialise the centroids to the eigengenes generated by wgcna, thus taking advantage of hc to carry out sensible initialization. implementation the standard wgcna procedure consider a gene expression profile matrix g n×m where n is the number of samples for a given condition, m is the number of transcripts and each g in g gives the quantification of the j-th transcript within the i-th sample. the standard wgcna procedure generates a squared adjacency matrix, between genes, based on their correlation. depending on whether the adjacency is signed or unsigned we will obtain networks either reflecting the direction of co-regulation or ignoring it, respectively. adjacency is defined as a d j=|c o r|β for genes i and j. the β parameter is an integer that modulates how smooth is the transition between the lowest to the highest possible co-regulation between genes. the wgcna methodology enables choosing β in such a way that the network shows a scale free topology property. this feature is commonly observed in biological networks. from the adjacency values, a new matrix with the same dimensions is created, the topological overlap matrix. this step alleviates the effect of noisy genes when obtaining the adjacency from correlation. once the network is built through the tom, it is converted to a distance matrix to use it as the basis for clustering. a dynamic tree-cutting algorithm is then applied to the dendrogram to generate a partition p={p …,p k} of disjunct sets of genes. thus, wgcna generates two main components which are useful for subsequent downstream analyses. on the one hand, the tom gives, for the j-th row/column, the level of co-expression of gene j with all of the genes in the network. the higher the value for a given pair, the tighter the interaction between them. furthermore, the sum of all row or column values for a gene, will give a measure of its overall level of co-regulation within the experimental condition, i.e. its ‘hubbiness’. thus, the tom is, in effect, the gcn. the other component produced by wgcna is the partition of gene sets, p, created from the tom. these partitions or modules often reflect cell types, common cellular functions or other biological subsystems reflecting, for example, immune function, or function related to the tissue under study. but the main utility of modules is to allow mapping gene groups to traits, when available. following the wgcna standard methodology, this is performed by looking for significant correlations between traits and the module ‘eigengene’. the eigengene summarizes the overall module activity in a given sample, and is obtained as the 1st pc component of the gene expression of genes belonging to the module. improvement of hierarchical clustering with k-means hc provides a convenient graphical representation of groupings that can be validated by biologists. one can readily obtain a suitable number of clusters from such an approach by ‘cutting’ the dendrogram at different heights, either manually or via various automatic algorithms. but, as we explained above, the final dendrogram strongly depends on how we measure distance between clusters. furthermore, once a gene falls under a subdendrogram, this decision cannot be modified under hc. if we consider how wgcna manages modules and eigengenes, it is assumed that each gene is highly correlated with other genes in its module. in other words, the module membership of the gene in its own module, measured as the pearson correlation between its expression and the module eigengene, should be higher than it is for any other module. however, we show here, from our real-data analyses, that 25% of genes would be better off in other modules. in this paper, we propose a post-processing step based on k-means to overcome all these limitations. it works on the partition p, leaving the tom unmodified. the k-means algorithm is well known and works on the n dimensional sample space of m points in an iterative fashion. it starts by setting a value for k, the number of clusters to discover and k centroids, one for each cluster. centroids are the representatives of each cluster, in such a way that a point g belongs to cluster i if the distance of such point to the cluster centroid is the minimum among all distances to all k cluster centroids. in standard k-means, given a partition of k modules, the the centroid for the i-th module c i={c i, …,c i,n} is generated as follows ci=1n∑j=1mgj,wheregi∈pi. however, in wgcna, the notion of a centroid is substituted by that of an eigengene. accordingly, our definition of k-means will use eigengenes as centroids. the concept of distance is a central element of k-means. it is important to note that distance in k-means is always defined between a point in the dataset and a centroid. euclidean distance is the the most commonly used distance in conventional k-means. however, given that we are constructing co-expression networks based on correlations, distance cannot be euclidean. modules should represent co-expressed genes instead. thus, and depending on the wgcna type of network, we should apply a distance between gene and eigengene based on the co-expression measure used. we will limit our discussion to signed networks. these specific types of networks will separate up- from down-regulated genes in different modules, which is usually of biological interest. they are also convenient for downstream analysis as correlation of genes and eigengenes will be positive, which eases a posteriori analyses. in signed networks, wgcna uses co=12), as a normalised measure of co-expression between the expression profile of a gene g i and a eigengene e g j, where by default c o r() is the pearson correlation coefficient. accordingly, we use 1−c o as distance. it is worth noting that hc needs a distance matrix between all genes, i.e. 1−t o m. k-means needs instead a computable distance definition between gene and eigengene. finally, on the basis of this definition of centroid and distance, genes are reassigned to the partitions induced by the new centroids, iteratively. if a stopping criterion is met, the algorithm finishes. otherwise, a new iteration is performed. we note that wgcna is computationally optimized to use pearson correlation. other correlation measures are in principle possible, including spearman’s rank correlation coefficient. however, in our own investigations we observe an increase in computation time of at least × when using spearman correlation without seeing any conclusive improvement with respect to the biology of networks. thus, throughout this paper we perform analyses using pearson correlation. we propose a general procedure which obtains, from a g n×m matrix of gene expression profiles from n samples and m genes, a clustering partition p of such genes by incorporating the standard wgcna process together with a post-processing of the partition obtained from it. the original contribution of this paper is described in steps from to below. step 1: initialization. let g n×m be a dataset of n samples and m genes for a given condition. let d() denote a distance function between a gene in g and an eigengene. let f c() denote a function which takes a clustering partition p={p …,p k} as an argument and generates centroids step 2: β=wgcna::picksoftthreshold$powerestimate step 3: obtain a tom, given g and β step 4: generate a partition p hc={p …,p k} with 1−t o m as a distance matrix and with average linkage hierarchical clustering and dynamic cutting height. step 5: let c={c …,c k} be a set of k vectors of n components which denote the centroids of the k-means clustering. step 6: initialize c with f c step 7: create a new partition p km={p …,p k} by assigning each gene g i, 1≤i≤m to a p j∈p km such that d≤d, 1≤t≤k holds. step 8: if the termination criterion holds, then stop. if not, generate new centroids c with f c and go to step note that, in this algorithm, we left d() and f c() undefined. however, within this paper we define d() according to eq. and f c() as the module eigengene. computational complexity of the proposed approach conventional wgcna gcn construction needs three sequential steps: obtaining the soft threshold to account for scale free topology, which has a computational complexity that depends on the number of genes and samples, obtaining the tom matrix, which has a complexity o, where n is the number of genes, as it has to construct a n×n squared matrix of adjacencies between all n genes, and hierarchical clustering, which has a complexity o. overall, wgcna’s computational complexity is o. wgnca’s space complexity is also o because it needs to maintain the tom in memory for hc to get the clusters. the computational complexity of k-means fits well with wgcna’s complexity. its time complexity is o where n is the number of genes, k is the number of clusters and it is the number of iterations. assuming that k,i t< using it as a post-process is very affordable in terms of computation time. note that k-means does not require the tom matrix in memory as the only distances it requires are between genes and eigengenes, and these we obtain on the fly by using eq. stopping criterion it is reasonable to assume that a sufficiently high number of k-means iterations will always be able to decrease the number of misplaced genes to on the other hand, the algorithm’s time complexity means that it is possible to run a single k-means in a conventional laptop in a matter of a few minutes. this means that we could, in principle, design a stopping criterion based on the minimun number of misplaced genes being set to however, we note that a situation may exist where the algorithm may fall into an infinite loop without reaching the desired state. thus, the stopping criterion we include in the software package km2gcn tries to reach the desired value for misplaced genes but always within a limited number of iterations. we did not observe the mentioned infinite loop situation in any of our experiments.  we wished to assess the ability of our method to define gene groups that genuinely reflect biological function. this is non-trivial for the following reasons. firstly, many genes are known to be pleiotropic, i.e. a single gene can affect many traits. transcription factors are a good example of this but there are many other examples. by creating non-overlapping partitions we deliberately ignore this fact and implicitly assume a model in which genes are highly specialized. secondly, we are limited by technology and sample availability from producing optimum estimates of gene expression profiles. we therefore lack of all the necessary information to build the best model. finally, if we wanted to evaluate the functional similarity of genes within a module, again we do not know all functions that all genes may play in any condition. notwithstanding these caveats, we explored various approaches to provide a comprehensive and varied assessment of the effectiveness of our k-means hybrid method. in “materials and methods for the gcns used for our evaluations” section we describe the datasets used in our investigations and the particular pipelines used to obtain the corresponding gncs. in “dynamics of k-means when working on 1-tom distance space” section we show our hybrid approach works. in “is k-means doing a proper job?” section we digress to note that k-means actually optimizes the sum of squares of within cluster distance. in “k-means improves the ‘eigengene’ as a tool for analysis” section we show that the proposed approach improves modules as a tool for mapping with traits. in “k-means improves module preservation” section we suggest that k-means improves cluster similarity between conditions. in “k-means detects more accurate partitions than wgcna in simulated data” section we compare the accuracy of k-means against wgcna on simulated data. in “k-means improves functional enrichments” section we show that k-means improves a module’s functional characterization through well-known databases such as the gene ontology. finally, in “k-means improves brain specific cell type marker enrichment” section we present results that suggest that gene markers for specific cell types show a better arrangement in partitions generated from k-means. materials and methods for the gcns used for our evaluations we evaluated gcns in two well-known datasets. the first is focused on brain tissue exclusively and it is based on affymetrix human exon v microarray expression profiles from brain tissues. this dataset is well suited for evaluating the k-means extension to wgcna because it is well known, it comprises different brain regions and gcn networks created with the standard wgcna method have been published. the procedure used to create the gcns is as follows. sample outliers were identified by visual inspection after clustering the samples using hierarchical clustering with euclidean distance as the distance measure. the majority of the identified outliers had low interarray correlation, which is defined as the pearson correlation coefficient of the expression levels for a given pair of transcripts using all available data available. after outlier removal, the same process was repeated to check for additional outliers. the gcn constructed was of signed type, with β= for all tissues. using these settings, the hc wgcna partition was created using transcripts passing quality control. once the partition was created, additional transcripts were assigned to modules based on their highest module membership. each partition was refined afterwards with k-means. the second dataset is gtex, which is one of most comprehensive human datasets currently available for multi-tissue transcriptomics. the gtex v gene expression dataset comprises samples unevenly distributed across post-morten human tissues. we created networks for tissues. in sequential steps, starting from rpkm values of gene-level quantification provided by gtex, we selected all tissues with more than samples. for each tissue we retained ensembl genes with r p k m> seen in more than 80% of the samples. this produced a variable set of genes for each tissue, with a minimum of for skeletal muscle tissue and a maximum of for testis. we applied batch, gender, age and rin as known covariates for data correction and to account for unknown covariates we applied sva axes. for each dataset of filtered rpkm values, we applied the sva r package using svaseq() to generate sva axes. for network construction we used the residuals obtained by regressing the rpkm expression values with the known and unknown covariates with a generalized linear model. to construct the networks, we applied the algorithm introduced in “improvement of hierarchical clustering with k-means” section. note that the differences between the ukbec and gtex networks are important and makes them well suited and complementary for the purpose of our study. the ukbec gene expression dataset is microarray based, while the gtex gene expression dataset is based on rna-seq technology, with rpkm quantification. ukbec networks are restricted to brain tissues while gtex networks cover tissues, including brain tissues. in summary, we have gtex gncs, brain specific ukbec gcns, gtex sample sizes range from n= to n=; ukbec sample sizes range from n= to n=; we have a variable number of genes used in the gtex gcns, in the range to for skeletal muscle and testis respectively; and the same probes for all ukbec gcns. finally, note there is a much higher variability in the number of modules per gcn in gtex, than in ukbec,. please note that throught the paper we use abbreviations to refer to tissues. please see table for the correspondence between abbreviations and brain region names. dynamics of k-means when working on 1-tom distance space as outlined in “improvement of hierarchical clustering with k-means” section, our proposed algorithm does not modify the distance matrix but acts later, on the partition p={p …,p k} taking k from the number of modules discovered by the hc used within wgcna. k-means acts iteratively creating centroids from modules and deciding for each gene, on the basis of the new centroids, which one is nearest to the gene. if, in the current iteration, the gene is nearest to a different centroid, then the k-means algorithm assigns it to the corresponding module. thus, in each iteration a new partition is generated with the changes applied to the former partition. i∈p, on average 30% of other modules within the gcn contribute with genes to its final gene set configuration. genes that leave their hc module have a module membership at that module of on average, with standard deviation of. genes arriving to a module for the final k-means partition show an average mm on arrival,, with standard deviation. fig. upper plot shows the evolution of the number of moved genes between any pair of modules p i and p j across k-means iterations for ukbec-microarray dataset. bottom plot shows the average module membership of genes moved across iterations for the ukbec-microarray dataset in comparison with average module membership for all the genes the lower panel of fig. focuses on the ukbec gcns, and on how mm evolves with iterations. dashed lines show, for each tissue, the average mm of moved genes, at each iteration, defining mm with reference to the original wgcna partition. initially, the algorithm focuses on moving genes with very low mm, but following this it focuses on genes with higher mm and then stabilizes. the solid lines show the average mm of all the genes in the network for each iteration. this dramatically increases over the first iterations and then smoothly and monotonically increases across additional iterations. this suggests that, over time, genes’ assignment to a module becomes stronger. is k-means doing a proper job? k-means is designed to optimize the sum of squares of within clusters distance. we define the within cluster distance, denoted with w, for a partition p as w=∑k=1k∑p=k||gi−ck|| where k is the number of clusters within p, p refers to the cluster that g i belongs to, and c k is the eigengene for the k-th cluster. we could, alternatively, define the distance between clusters, b. note that for a given set of genes, we can obtain the sum of distances between all gene pairs. if we denote this measure by t for a given gene expression profile g, we can decompose it into t=b+w for any given p obtained from g. this means that maximizing the between-clusters distance is equivalent to minimizing the within-cluster distance. we observe that k-means monotonically decreases w increases) across iterations. the algorithm generates a higher w at the early iterations, which decreases to a lower level in later iterations. this behaviour is in line with the shape of the gene changes curves of the upper plot in fig. higher number of moved genes imply higher decreasing rate of within cluster distance. fig. the within cluster distance evolution during the k-means runs for the ukbec datasets this behaviour is also in accordance with what we see in fig. this plot shows, for each module p i∈p, and the specific case of ukbec’s cerebellum cgn, the distance between the eigengenes for the same module, as they are created during successive iterations. over time, the eigengene vectors stabilize across iterations suggesting that cluster definition becomes stable. fig. euclidean distance of successive module eigengenes along the k-means iterations for cerebellum samples for ukbec datasets k-means improves the ‘eigengene’ as a tool for analysis one of the main applications of wgcna partitions is searching for associations between gene clusters and traits. traits are usually given as a vector of n components where n is the number of samples. on the other hand, modules are comprised of m vectors of the same form where each vector is the expression profile for the corresponding gene. to assess the correlation between the trait and the gene module, wgcna transforms module p i into an eigengene. from this, the correlation between the trait and the eigengene can be easily obtained. there are several applications for the eigengene. for example, it can be used to provide a measure of how strong is the membership of each gene g∈p i to the i-th module, by correlating its expression with the eigengene, resulting in the mm of g. let this module membership be denoted with m for gene g and module i. it is assumed that a good p would be one such that, the number of genes g with m<m when g∈p i, for any i≠j, is low. let us call such genes ‘misplaced genes’. we would prefer a partition in which the number of misplaced genes is minimum. to assess the number of misplaced genes after applying k-means, we performed investigations in the ukbec and gtex tissues. in ukbec, using k-means with only iterations, we get a maximum of misplaced genes in the putamen and a minimum of for occipital cortex and a mean of misplaced genes per partition. in the wgcna partitions, the maximum number of misplaced genes is in temporal cortex, the minimum is in white matter and an average of genes; times more than with the k-means algorithm. in gtex clustering partition modules, the average number of misplaced genes in modules from a wgcna partition is after applying k-means, it is only. k-means improves module preservation one component of wgcna provides a convenient tool for the analysis of module preservation. given a partition, p, constructed from a network obtained from a given set of samples s, we can test whether the features of each module p i∈p are preserved in an alternative set of samples s ′. preservation analysis is based on estimating, for some statistic of interest, differences between what is observed and what is obtained by random permutation. for example, one statistic of interest is the gene correlation with the eigengene. through a simple transformation one can check whether the values obtained in the reference network are maintained for the same genes within the other network. wgcna uses the ‘z-summary’ statistic as a general summary of all the different statistics used. to assess the effect of k-means on z-summary, we performed the same investigations on both ukbec brain tissues and on the gtex brain tissues. note we focus on brain tissues within gtex as comparison of preservation only makes sense for tissues that are similar. within each ukbec and gtex tissue gcn, we compared the preservation of all the partitions generated by wgcna alone with the preservation obtained by applying k-means to each of them. a permutation analysis on tissues generates, for each tissue t, and for each module p i within the corresponding tissue network, a vector of z-summary values corresponding to the preservation of p i in the other tissues. table number of new modules from a tissue that are preserved on another tissue after applying the k-means to the standard wgcna partitions from table it is apparent that there is an overall increase in the number of modules preserved under k-means. in the ukbec gncs, there is an improvement in cases, no improvement in cases, and only case with a worse preservation. the average improvement in modules preserved for ukbec is. in the gtex gcns, there is an improvement in module preservation in cases, no improvement in only cases and a decreased preservation in just cases. the average number of modules improved by the k-means method is. this suggests that k-means creates less noisy modules as similarities between tissues are more apparent. finally, it is worth noting that each tissue is expected to have specific modules, i.e. modules that will be poorly preserved in other tissues because they are exclusive from that tissue, reflecting study-specific or sample-specific gene subsystems. k-means detects more accurate partitions than wgcna in simulated data we wanted to test whether k-means improves the accuracy of partitions p with respect to those obtained under standard wgcna. to this end we investigated networks based on ‘synthetic data’. the wgcna package provides a gene expression profiling simulation method simulatedatexp(), which is a convenient method for generating artificial data sets that mimic the properties of real datasets. the simulation method works with the eigengene of gene expression for each gene belonging to the module. the simulation method requires, as arguments: a matrix with the eigengene for each module; the proportion of the total gene pool that one will find within each module; and the number of genes to be simulated. note that the number of samples we want to simulate appears implicitly as the length of each eigengene. the method returns two elements: a gene expression profiling data set, let us denote it with d, that we can use to construct gcns; and the ideal clustering partition of the simulated gene expression profiling, here denoted by p. thus, if we rely on the effectiveness of this simulation method, then a simulated data set d, we will prefer a gcn construction algorithm a to algorithm b if the distance between a and p is smaller than between b and p, where a and b are the clustering partitions we get after constructing gcns on d with a and b, respectively. the accuracy of an algorithm a is defined by the similarity of the theoretical optimal partition within the synthetic data to the partition constructed by a. in order to test whether k-means performs any better than standard wgcna on simulated data, we constructed a plausible set of simulated gene expression profiles. we used gtex and test with them both k-means and standard wgcna on gcn construction. the accuracy of an algorithm a will be defined as how similar are the theoretical optimal partition within the synthetic data, and the partition constructed by a. in order to test whether k-means performs any better than standard wgcna on simulated data, we constructed a plausible set of simulated gene expression profiles. we used gtex standard wgcna gcns as the simulation seed for the generation of a synthetic gene expression profile. we focused on the gtex dataset rather than ukbec, because the gcns comprise a usefully heterogeneous network dataset. the simulated data process produced a gene expression profile and a theoretical ideal clustering partition for such profile. we used this theoretical ideal partition to evaluate standard wgcna and k-means accuracy. to estimate accuracy we use three different statistics: the rand index, also implemented within wgcna, the jaccard coefficient and the similarity index, all of them implemented within clv r package. results for all the experiments appear in additional file each row corresponds to a gtex tissue, the randsimvswgcna column corresponds to the rand index between the ideal partition and that obtained with wgcna on the simulated data. the randsimvskm column corresponds to the same index when using k-means. the other four columns correspond to the jaccard coefficient and the similarity index. the k-means refinement generate higher values in all the cases for all three indexes. these results are illustrated in fig. fig. results on performance of standard wgcna and k-means on simulated data sets that used the gtex wgcna gncs as seed for simulation. we display the same results using three different indexes of similarity between cluster partitions. the k-means method outperforms standard wgcna with all three indexes used k-means improves functional enrichments the gene ontology is a curated database for gene annotation which can be used for the functional characterization of gene sets. given a set of genes, and a subset of those genes, an enrichment analysis can be performed on go annotations to search for terms in the ontology that are significantly enriched in the subset of genes relative to the full set. the number and strength of significant terms obtained in this way can be used to measure the biological functionality of the module. given two different partitions p and p ′ created from the same tom, we would prefer the partition that generates more significant go terms if we assume go to reflect a biological ground truth as this would suggest that the preferred partition makes more biological sense. we used the gprofiler r package to obtain enrichment p-values, avoiding eia terms in go and requiring a correction for multiple testing with gscs, as developed by the authors of the package. we describe below a series of investigations to characterize the improvement in a module’s biological functionality. global annotation term significance consider a partition p={p …,p k} of genes arranged into modules p i, 1≤i≤k. now suppose we want to perform a gene set enrichment analysis on each p i∈p based on the gene ontology. go is a list of ontological terms, organised into three main branches: bp, mf and cc. genes within the database will be associated with a number of terms from each branch. thus, for each term in go, and given the list of genes in p i, we can apply a contingency test, e.g. fisher exact test, under the null hypothesis that the genes in p i show no significant overlap with the set of genes associated with the term. with an appropriate correction for multiple testing, we define as significant the association of the list of genes in p i with the corresponding term, when the corrected p-value is <. we then aggregate all these p-values for a module in a single measure of significance as follows. for each p i∈p, we use sgo=∑pvaluej∈test−log where t e s t is the set of p-values, p v a l u e j, of significant terms associated with the genes in partition p i, emerging from the analysis. in this way, sgo=∑pi∈psgo can be used to aggregate all the biological signals of a whole partition p. given a choice of partitions, we prefer p to p ′ when s go>s go. a displays for each ukbec and gtex gcn, the relative improvement between the standard wgcna partition, p and the k-means partition, p ′ by sgosgo− fig. the left plot’s light blue blue bars show the percentage of relative improvement by k-means with respect to wgcna s go statistic. values in red are those that k-means fails to improve. the right plot shows cell type enrichment improvement in the same way, for the ukbec gcns and the gtex brain networks. again, values in red are those that k-means fails to improve the average improvement is 13%. overall, there is improvement in all ukbec tissues and in out of gtex gcns, and the overall improvement is significant. does a higher enrichment implies less informative modules? high values of the s go index are of interest, as we prefer a partition p over p ′ if s go>s go. however, it is possible that modules show better s go values after k-means because the module have more annotation terms that are generic, and therefore less descriptive about the specifics of the tissue studied. in order to assess this, we applied the notion of information content. we used the gosim package which applies information-based metrics to gene ontology terms. the metric i c for a term t belonging to an ontology is defined as: ic=−logp, where p is the probability of observing t within the annotations available within that ontology. ideally, we prefer modules with more go terms, which are more significant and more informative. from previous sections we know we have more significant networks thanks to s go. but is k-means capable of not only improving significance but also of maintaining the level of information of the modules if not increasing it? fig. relation between frequency of appearance of go annotation terms across all gtex gcns and ic. terms appearing more times tend to have lower ic. regression lines show that k-means gets better ic values for highly repetitive terms is the overall ic obtained by k-means degraded as a consequence of obtaining more significant terms per gcn in comparison with standard wgcna? to assess this we regressed the information content of the significant annotation terms against the frequency of appearance in the gcn annotation sets. we found a tendency towards higher ic in k-means gcns. this suggests that k-means annotations are more specific, and therefore more useful. is the increase in enrichment better than random? in “dynamics of k-means when working on 1-tom distance space” section we noted that one of the changes within cluster partitions after applying k-means is that module sizes change and many modules will increase their size considerably. it is fair to assume that modules increasing their size in genes, will also increase their s go enrichment. there is a significant pearson correlation between increase in module size after k-means and increase in number of significant annotation terms. the question arises, therefore, of what is the real contribution of k-means in comparison to a random shift of genes between modules? in order to answer this question, for each of the gtex tissues and their corresponding wgcna and k-means partitions, we identified those genes that were changed at the wgcna partition to create the k-means one. then, in a single step, we randomly assigned these genes to other modules in such a way that we kept the same module sizes obtained with k-means. via this algorithm, we produced new partitions in which the genes that remained unchanged from wgcna to k-means stayed at the same modules, but those genes that were changed by k-means were again changed but this time in a random fashion. go statistic. plot shows the number of significant go terms. fig. effect of random assignment of genes selected by k-means, on a wgcna partition, to be changed from one module to another. plot refers to s go values and to number of significant terms in 89% of the modules, k-means finds the same number or more significant go terms than the random placement of misplaced genes. 88% of the final modules show equal or better s go index using k-means compared to random. aggregating the results by tissue, k-means placement performs better in all the cases. interestingly, the random placement of misplaced genes prevents enrichment at the wgcna partition it starts with. this is important because even though many genes at the wgcna partition are not touched by the random approach, moving genes randomly will nevertheless worsen these genes’ functional annotations. this suggests that both the number of significant terms and the s go index have a reasonable sensitivity. k-means improves brain specific cell type marker enrichment one interesting property of wgcna gcns is that partitions created from them can be useful when studying cell-type specific gene networks. in studies where samples come from bulk tissue, it is most likely that these samples will be comprised of different cell material. in consequence, the gene expression profiles obtained from them should reflect this heterogeneity in some way. wgcna’s gcns handle this heterogeneity in an elegant and convenient manner: they often generate gene clusters within partitions which are specialized on a given cell type, i.e. they present a highly significant enrichment of markers for a given cell type. we wanted to assess the effects of k-means on this particular feature. to do this, we used three different resources defining cell-type specific gene sets. these were wgcna’s brain lists,, and two alternative brain specific sources, labelled here external and cahoy. we evaluated each partition’s modules from the ukbec gcns and the brain tissue gtex gcns, using both standard wgcna and k-means. this evaluation generated two matrices of p-values, with each gene dataset in a row and each specific module from each of the networks in a column. p-values reflect a fisher’s exact test for whether there is significant concentration of the corresponding gene sets in the tested module. we include in additional files and the results for standard wgcna and k-means, on the ukbec gcns. note that, in these plots, both columns and rows have been clustered based on −l o g so it can be better seen how modules from different tissues cluster together at columns, and also how different gene sets cluster among rows. these heat-maps reveal strongly clustered areas corresponding to groups of cell-type specific genes sets within most, if not all, of the tissues. more specifically, we see four groups cell-type specific gene sets corresponding to microglia, astrocytes, oligodendrocytes and neurons. in the ukbec k-means heat-map, using a significance cut-off of 10−, almost 65% of the modules show cell type enrichment. within these, modules show a single cell type signal. in the wgcna heat-map, 63% of the modules show cell type enrichment, with showing single cell type signals. b compares the two enrichment matrices, by aggregating all the cell-type enrichment −l o g transformed p-values as we did for the gene ontology enrichment in “k-means improves functional enrichments” section. each bar represents the sum of all values of the corresponding heat-map, for modules of the given tissue. according to this, we always see an improvement in ukbec networks and in most of the the gtex networks. the overall improvement is significant.  our study shows that an additional k-means step, when used as an adjunct to wgcna, improves the partitions generated from gene co-expression networks. our method is not an alternative to wgcna, instead it is an additional step to the standard wgcna pipeline. indeed, our method can be applied to any general hierarchical clustering algorithm, and as such it could be usefully applied to any hierarchical clustering based approach for network generation, not just gene co-expression networks. we evaluated our method using two contrasting gene expression datasets representing a variety of different tissues, one obtained with microarray technology, and the other with rna-seq. using a variety of approaches, we demonstrate improved performance of our k-means method in both datasets. furthermore, we also demonstrate improved performance using simulated data generated from the gtex dataset. we show via these analyses that it is possible to obtain better partitions for the same networks via our k-means method. our method generates modules with fewer misplaced genes with respect to their eigengene, and this implies that the eigengene is a better representative of the phenomena hidden behind the particular set of genes belonging to the module. using gene ontology enrichment analyses, we also show that our partitions are enriched for biological functionality. statistically significant s go enrichment is seen in all ukbec cgns and in out of the gtex gcns. our partitions have improved modules preservation, which also suggests that the clustering is more accurate from a biological point of view. although some gene modules are specific of each tissue, it is a desirable property of most gcn partitions to be highly replicable under the assumption that a preserved module is more likely to be a genuine module. our analyses suggest that k-means favours the creation of more genuine modules and these results are seen in both ukbec and gtex gcns. our k-means method also creates partitions in which gene sets representative of specific brain cell types are seen in modules with increased statistical significance. this suggests, once again, more biologically genuine modules. gcn construction is likely to become an increasingly important analysis, as genomics and transcriptomics are increasingly applied to aid clinical diagnosis and prognosis. methods that generate more reliable and robust gene gcns will enable improved prediction of inter-gene relationships and gene function, with a variety of applications. availability and requirements ukbec data has accession code gse all information about tissues, samples and quality control can be found there. gtex rpkm gene expression v was used in this paper and downloaded from the gtex portal: http://gtexportal.org/home. regarding the software we present here, this is the availability and requirements. project name: km2gcn project homepage: https://github.com/juanbot/km2gcn operating system: linux/windows/mac programming language: r other requirements: wgcna r package and gprofiler license: lgpl additional files additional file lists tissues, samples and genes used for the creation of each gcn. additional file includes comparative results of standard wgcna and k-means on simulated data. additional file the second one corresponds to k-means. both on ukbec datasets. values higher than are set to colors at the top of columns correspond to tissues, the tissue legend is at the bottom of columns and cell marker gene set used on the right side. additional file heat-maps showing −l o g from ficher’s exact test on significant concentration of specific cell marker gene sets on each tissue module. the one within the additional file corresponds to the standard wgcna. abbreviations bpbiological process category of the gene ontology cccellular component category of the gene ontology eiaelectronic inferred annotations gbaguilty by association gcngene co-expression network gogene ontology gtexgenotype-tissue expression project hchierarchical clustering mfmolecular function category of the gene ontology mmmodule membership of a gene sftscale free topology tomtopological overlap matrix ukbecunited kingdom brain expression consortium wgcnaweighted gene-co expression network 
weighted gene co-expression network analysis is a widely used r software package for the generation of gene co-expression networks. gene co-expression networks on braink-means applied to wgcnaassessment of better gene clusters on bulk tissuehttp://dx.doi.org/ /501100000265medical research councilmr/k01417x/1hardy john alzheimer’s research uk issue-copyright-statement© the author 2017 we observed substantially improved module properties, including: few or zero misplaced genes; increased counts of replicable clusters in alternate tissues; improved enrichment of gene ontology terms improved cell type enrichment signals; and more accurate partitions in simulated data according to a range of similarity indices. the results obtained from our investigations indicate that our k-means method, applied as an adjunct to standard wgcna, results in better network partitions. we propose k-means clustering as an additional processing step to conventional wgcna, which we have implemented in the r package km2gcn. wgcna generates both a gcn and a derived partitioning of clusters of genes. these improved partitions enable more fruitful downstream analyses, as gene modules are more biologically meaningful.
 weighted gene co-expression network analysis is a widely used r software package for the generation of gene co-expression networks. wgcna generates both a gcn and a derived partitioning of clusters of genes. we propose k-means clustering as an additional processing step to conventional wgcna, which we have implemented in the r package km2gcn.  we assessed our method on networks created from ukbec data, on networks created from gtex data, and on simulated networks derived from gtex data. we observed substantially improved module properties, including: few or zero misplaced genes; increased counts of replicable clusters in alternate tissues; improved enrichment of gene ontology terms improved cell type enrichment signals; and more accurate partitions in simulated data according to a range of similarity indices.  the results obtained from our investigations indicate that our k-means method, applied as an adjunct to standard wgcna, results in better network partitions. these improved partitions enable more fruitful downstream analyses, as gene modules are more biologically meaningful. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords gene co-expression networks on braink-means applied to wgcnaassessment of better gene clusters on bulk tissuehttp://dx.doi.org/ /501100000265medical research councilmr/k01417x/1hardy john alzheimer’s research uk issue-copyright-statement© the author 2017 