Gene co-expression analysis for functional classification and gene–disease predictions

1A key objective in biological research is to systematically identify all molecules within a living cell and how they interact. However, the functions of many genes are still not understood, a situation that has only become more complex with the recent identification of many novel non-coding genes. With the development of high-throughput technologies including microarrays and RNA sequencing (RNA-seq), and their respective data-analysis methods, the functional status of a gene can now be identified from a systematic perspective. One method to infer gene function and gene–disease associations from genome-wide gene expression is co-expression network analysis (Figure 1), an approach that constructs networks of genes with a tendency to co-activate across a group of samples and subsequently interrogates and analyses this network.

2Gene co-expression networks can be used for various purposes, including candidate disease gene prioritization, functional gene annotation (Figure 1) and the identification of regulatory genes. However, co-expression networks are effectively only able to identify correlations; they indicate which genes are active simultaneously, which often indicates they are active in the same biological processes, but do not normally confer information about causality or distinguish between regulatory and regulated genes. An increasingly used method that goes beyond traditional co-expression networks is differential co-expression analysis. This approach identifies genes with varying co-expression partners under different conditions, such as disease states, tissue types and developmental stages, because these genes are more likely to be regulators that underlie phenotypic differences. The regulatory roles of such genes can be further investigated by integrating data types such as protein–protein interactions, methylome data, interactions between transcription factors (TFs) and their targets, and with sequence motif analysis of co-expressed genes. This aids in the identification of regulatory elements such as TFs, expression quantitative trait loci (eQTLs) and methylation patterns that affect the expression and composition of co-expression modules.

3Gene expression and regulation can be highly tissue-specific, and most disease-related genes have tissue-specific expression abnormalities. The increased availability of expression data for multiple tissues has allowed for differential co-expression analysis, which can identify both tissue-specific signatures and shared co-expression signatures. These tissue-specific signatures can be disrupted in tissue-specific diseases and would not be detected in analyses aggregating multiple tissues. Even when no sample classification is available, subpopulation-specific modules can be resolved, an approach that has been particularly successful in classifying different cancer subtypes to provide prognostic markers. Differential co-expression analysis is also useful for analysing data sets in which the subpopulations are unknown, e.g. large-scale single-cell RNA-seq data. While differential co-expression methods are sensitive to noise, they are becoming more effective with the increase in RNA-seq data quantity and quality. RNA-seq further permits co-expression analysis to focus on splice variants and non-coding RNAs.

4In this review, we provide an introduction and overview of what constitutes a co-expression network, followed by a guide of the different steps in co-expression analysis using RNA-seq data. We then describe commonly used and newly emerging methods and tools for co-expression analysis, with a focus on differential co-expression analysis to identify regulatory genes that underlie disease. We conclude with a discussion of the integration of co-expression networks with other types of data, to e.g. infer regulatory processes, and with future prospects and remaining challenges in the field.

5A co-expression network identifies which genes have a tendency to show a coordinated expression pattern across a group of samples. This co-expression network can be represented as a gene–gene similarity matrix, which can be used in downstream analyses (Figure 1). Canonical co-expression network construction and analyses can be described with the following three steps.

6In the first step, individual relationships between genes are defined based on correlation measures or mutual information between each pair of genes. These relationships describe the similarity between expression patterns of the gene pair across all the samples. Different measures of correlation have been used to construct networks, including Pearson’s or Spearman’s correlations. Alternatively, least absolute error regression or a Bayesian approach can be used to construct a co-expression network. The latter two have the added benefit that they can be used to identify causal links and have been explained elsewhere. For a discussion of other types of similarity measures, we refer to. Many of these similarity metrics can also be used to construct protein–protein interaction networks, which were compared using cancer data in.

7In the second step, co-expression associations are used to construct a network where each node represents a gene and each edge represents the presence and the strength of the co-expression relationship (Figure 1).

8In the third step, modules (groups of co-expressed genes) are identified using one of several available clustering techniques. Clustering in co-expression analyses is used to group genes with similar expression patterns across multiple samples to produce groups of co-expressed genes rather than only pairs. The clustering method needs to be chosen with consideration because it can greatly influence the outcome and meaning of the analysis. Many clustering methods are available, including k-means clustering and hierarchal clustering, and are discussed in detail in. Modules can subsequently be interpreted by functional enrichment analysis, a method to identify and rank overrepresented functional categories in a list of genes.

9In co-expression analysis, it is important to consider the heterogeneity of the samples. Tissue-specific or condition-specific co-expression modules may not be detectable in a co-expression network constructed from multiple tissues or conditions because the correlation signal of the tissue/condition-specific modules is diluted by a lack of correlation in other tissues/conditions. However, limiting co-expression analysis to a specific tissue or condition also reduces sample size, thereby also decreasing the statistical power to detect shared co-expression modules. Therefore, methods that do not distinguish between tissues or conditions should be used for identification of common co-expression modules, while differential co-expression comparing different conditions or tissues will be better for identifying modules unique to a specific condition or tissue.

10In a correlation-based co-expression network, correlation measures have values between −1 (perfect negative correlation) and 1 (perfect positive correlation). In an unsigned network, the absolute correlation values are used, which means that two negatively correlated genes will be considered as co-expressed. This causes negatively correlated genes to group together. Because those genes are likely to be also positively co-expressed with a completely different set of genes, these genes also group into the same module and disrupt the structure of the network. A signed network solves this problem by scaling the correlation values between 0 and 1 so that values <0.5 indicate negative correlation and values >0.5 indicate positive correlation. A signed method creates networks where biologically meaningful modules (such as those representing a specific biological process) are better separated. Thus, a scaled value close to 0 indicates negative correlation, a feature which may be particularly interesting when microRNAs (miRNAs) are incorporated into the network, as these are known to exert their function mainly through down-regulation of other genes. This also holds true for some long intergenic non-coding RNAs (lincRNAs).

11In a weighted network, all genes are connected to each other, and these connections have continuous weight values between 0 and 1 that indicate the strength of co-regulation between the genes. In an un-weighted network, the interaction between gene pairs is binary, i.e. either 0 or 1, and genes are either connected or unconnected. An un-weighted network can be created from a weighted network by, for example, considering all genes with a correlation above a certain threshold to be connected and all others unconnected. We focus on weighted networks in this review because (to date) they have produced more robust results than un-weighted networks.

12Co-expression networks can be constructed from gene expression data obtained from microarray or RNA-seq technology. One of the major benefits of RNA-seq is that it quantifies the expression of the over 70 000 non-coding RNAs not usually measured with microarrays, including recently annotated lincRNAs, many of which are thought to have regulatory roles and to play a role in disease. Therefore, to gain a better understanding of the regulatory mechanisms driving biological processes, non-coding RNAs need to be considered in analyses.

13RNA-seq also has other benefits. It increases accuracy for low-abundance transcripts, has a higher resolution for identifying tissue-specific expression and distinguishes expression profiles of closely related paralogues better than microarray-derived profiles. RNA-seq can also distinguish between the expression of different splice variants, which can have distinct interaction partners and biological functions. Co-expression analysis on RNA-seq data can assign putative roles to these splice variants and lincRNAs, and identify diseases in which they might play a part. A limitation of co-expression analysis on the splice variant level is the introduction of biases because it is difficult to determine which splice variant is expressed if multiple splice variants share the same expressed exon.

14As an example of RNA-seq's utility with isoform- and exon-specific expression level measurements, exon-level expression was used to construct a co-splicing network. In a gene co-expression network, expression of different transcripts originating from the same gene is usually aggregated, which can lead to biased co-expression signals. In a co-splicing network, this issue is resolved by considering the exon-expression-level distributions within a gene when calculating gene co-expression correlation. In biological terms this means that the expression of two genes is only considered to be correlated if their different splice variants show co-ordinated expression. If this is not the case, they are not considered to be co-expressed even if the overall expression levels of the genes are correlated. This approach has identified novel functional modules, which would not be detected using traditional co-expression networks. Additionally, genes that contain multiple exons and transcripts acquired more relevant positions in the network using this method, a reassuring result given that splice variants can have different functions and are thus likely co-expressed with functionally distinct partners, which co-splicing networks account for.

15A different approach is to determine the expression of different isoforms originating from the same gene based on the distributions of reads mapping to its various exons. This method is used by SpliceNet, which effectively divides the reads mapping to an exon shared with two isoforms proportionally to the total expression of each of the two whole isoforms. This means that if two isoforms, isoform A and isoform B, share only one exon X (to which a number of reads map), but there are no reads mapping to the other exons of isoform A, whereas some reads map to the exons of isoform B, all reads mapping to exon X are then assigned to isoform B, resulting in isoform A being considered as not expressed at all. Although this elegant solution was validated using simulations, no experimental validation was conducted.

16The most common way of constructing RNA-seq-based co-expression networks is to merge all overlapping gene isoforms in the RNA-seq data analysis and then construct the network at the gene level. This approach, however, loses information about different transcripts encoded by the same gene. Alternatively, transcript-based co-expression networks can be constructed. The drawback of these networks is their dramatic increase in size owing to the many gene isoforms and non-coding RNAs. As co-expression networks are square matrixes, the size of the network increases quadratically (n2) with the number of genes included. As there are ∼200 000 annotated transcripts in the human genome (according to Ensembl GRCh38.p5 (human) annotation) and only ∼20 000 protein-coding genes, the resulting network increases 100-fold in size, greatly increasing the computational resources needed for the analysis. One solution to this problem is to build co-expression network blocks from subsets of the data and combine these blocks at a later point in the analysis. We recommend users to be cautious with block-wise clustering, however, as it may influence the results of subsequent module detection analyses, and it is unclear how well these perform when large numbers of blocks are used.

17RNA-seq analysis entails multiple steps that include obtaining expression estimates from the sequenced reads, data normalization and quality control. Different tools and methods to obtain reliable expression counts from RNA-seq data were recently reviewed in, and these will not be reviewed here.

18In our experience, different normalization methods introduce different biases in co-expression analysis, usually towards positive correlation. New methods are continuously being created to tackle these normalization issues. The recently published method extracting patterns and identifying co-expressed genes (EPIG) from RNA-Seq data (EPIG-seq), for example, is designed to calculate gene correlation across RNA-seq samples, being unaffected by read-depth differences between samples and the large abundance of 0 values present in RNA-seq-derived expression matrices. Biases originating from the large abundance of 0 values are even more pronounced in single-cell experiments because of low RNA quantities per cell. Specific tools have been created for analysis of single-cell RNA-seq data and are reviewed in. Although some studies comparing different normalization methods for RNA-seq data are available, more comprehensive comparison studies incorporating newer methods are needed.

19To create co-expression networks from RNA-seq data, a 20-sample minimum has been suggested, and increased sample sizes produce networks with a higher functional connectivity. Not surprisingly, higher quality data tend to result in more accurate co-expression networks. It is therefore essential to set cut-off thresholds for data quality control. A higher total read depth for RNA-seq samples increases the accuracy of the expression measurements, especially for genes with low expression. For RNA-seq data, sequencing depth cut-off thresholds are usually selected arbitrarily. Several co-expression studies have used a cut-off of 10 million reads per sample. Co-expression networks constructed using this cut-off have been suggested to have a similar quality to microarray-based co-expression networks if constructed from the same number of samples, but decreasing in quality with fewer reads. The percentage of mapped reads is another frequently considered cut-off in which samples with <70% or 80% of the reads mapping to the genome are removed. Giorgi et al. demonstrated, using 65 Arabidopsis thaliana samples with 12 million reads but applying only a 30% mapping cut-off threshold, that the resulting RNA-seq-based co-expression network had a lower similarity to biological networks than microarray networks. Cut-off thresholds may vary per species, based on, among other factors, the quality of the genome annotation. As more and higher quality data become available, higher cut-off thresholds may be preferable.

20To ensure that a network is robust, bootstrapping can be used. This is the repetitive construction of networks by using random sets of samples (one sample can be part of multiple subsets) from the data, which are subsequently used to assess the reproducibility of the network created from the entire data set. Randomizing the data set (e.g. by randomly reassigning expression values to their gene/transcript identifiers and reconstructing the network) can also help identify correlations that occur stochastically because of specific biases rather than as a result of biologically relevant interactions.

21Clustering is used to group genes that have a similar expression pattern in multiple samples. The resulting modules often represent biological processes and can be phenotype specific.

22The most widely used clustering package for co-expression analysis is Weighted Gene Correlation Network Analysis (WGCNA). This easy-to-use tool constructs co-expression modules using hierarchical clustering on a correlation network created from expression data. Hierarchical clustering iteratively divides each cluster into sub-clusters to create a tree with branches representing co-expression modules. Modules are then defined by cutting the branches at a certain height (Figure 1).

23WGCNA was the first co-expression tool to be applied to RNA-seq data; it has effectively identified biologically relevant associations between phenotypes and modules, performing similarly to microarray-based analyses. An RNA-seq-based co-expression study on normal and failing murine hearts found that many lincRNAs are present in clusters correlating with the failing murine heart phenotype, suggesting a possible role of these non-coding RNAs in this disease. Co-expression analysis of RNA-seq data of lean and obese porcines identified obesity-related modules, and a link was found between obesity, the immune system and bone remodelling, with the study identifying CCR1, MSR1 and SPI1 as possible regulators in these processes. WGCNA was also used to identify biologically relevant associations from single-cell RNA-seq data. Regulatory mechanisms and genes underlying pre-implantation processes conserved between humans and mice were identified by using preservation detection defined by WGCNA, a feature that was later added to this package. Co-expression modules were identified for different developmental stages of human and mice separately. The modules identified for each stage were then compared between humans and mice to reveal a strong overlap between co-expression modules in oocyte formation in mice and oocyte and single-cell stage co-expression modules in humans. This suggests that humans and mice share core transcriptional programmes in early development, but diverge at a later stage.

24Co-expression modules identified by clustering are often large, and so, it is important to identify which gene(s) in each module best explains its behaviour. A widely used approach is to identify highly connected genes in a co-expression network (hub genes). Hubs are frequently more relevant to the functionality of networks than other nodes. This is also the case in biological networks, although mathematical derivations show that this is only the case for intra-modular hub genes (as opposed to inter-modular hub genes). Intra-modular hubs are central to specific modules in the network, while inter-modular hubs are central to the entire network (Figure 2). To identify hub genes, centrality measures, mainly ‘betweenness centrality’, are often used. Genes with high betweenness centrality are important as shortest-path connectors through a network. Connectivity is often used to measure network robustness and indicates how many genes need to be removed from the network before the remaining genes are disconnected. Identifying hub genes in co-expression networks has led to the identification of several genes essential in cancer, type 2 diabetes, chronic fatigue, other diseases and tissue regeneration.

25As there are usually multiple hubs or differentially connected genes in a module, it is not always clear which is the most important gene underlying a phenotype. Nor is there a guarantee that any of the hub genes is causal to a phenotype. One study using myocardial data from 1617 samples found that known foetal gene markers upregulated in gene networks common to developing and diseased myocardium were not hub genes. Another co-expression study in Salmonella found that hub genes are dispensable for growth, stress adaptation and virulence, suggesting that hub genes are not necessarily essential.

26A widely used approach to attach biological meaning to modules is to determine functional enrichment among the genes within a module using e.g. the tools described in Table 1. Assuming that co-expressed genes are functionally related, enriched functions can be assigned to poorly annotated genes within the same co-expression module, an approach commonly referred to as ‘guilt by association’ (GBA). GBA approaches are also widely used to identify new potential disease genes if a substantial proportion of the genes within a module are associated with a particular disease (Figure 1).

27When using a GBA approach it is important to remember that not every gene in a module necessarily correlates with a function or disease association for which it is enriched. Because co-expression modules often consist of a large number of genes, any overrepresentation of a functional process or group of disease-associated genes quickly becomes statistically significant, as often indicated by deceivingly low p-values. Misinterpretation of these low p-values may lead to the incorrect conclusion that all genes in a module play an important part in a particular process or disease. In reality, the fraction of genes in a module that relate to its main biological function is often <20%, and module-trait correlations can be relatively low (correlation < 0.5) even when statistically significant.

28Although there is ample evidence that co-expression analysis can help identify genes that play an important role in disease and biological functions, it remains difficult to infer causality from co-expression networks. Tools such as ARACNE and GENIE3 attempt to construct regulatory networks from co-expression networks. ARACNE removes indirect connections between genes (i.e. partners of a gene that have a stronger correlation with each other than with the gene itself), leaving only those connections that are expected to be regulatory. GENIE3 incorporates TF information to construct a regulatory network by determining the TF expression pattern that best explains the expression of each of their target genes. A limitation of GENIE3 is that TF information is required for it to perform better than random chance. The performance of these methods has been compared with gold standards defined by regulatory interactions experimentally validated in >150 studies. The comparison suggests that methods attempting to derive regulatory networks from co-expression networks alone can only reliably distinguish between true- and false-positive regulatory interactions if perturbation experiment data are used for network construction. A comparison between these tools and others, including WGCNA, showed that WGCNA and ARACNE perform best at defining the network structure of Escherichia coli, for which a well-defined regulatory network was used as a gold standard.

29Differential co-expression analysis can identify biologically important differential co-expression modules that would not be detected using regular co-expression or differential expression analyses. Genes that are differentially co-expressed between different sample groups are more likely to be regulators, and are therefore likely to explain differences between phenotypes. Differential co-expression analysis has been used to identify genes underlying differences between healthy and disease samples or between different tissues, cell types or species. Below, we provide an overview of commonly used and newly emerging methods and tools, separated into two categories: (1) approaches that identify differential co-expression between predefined sample groups (such as conditions, time points or tissue types) and (2) approaches that do not require prior knowledge about sample groups and use an algorithm that identifies co-expression clusters in a priori unknown subpopulations of the samples.

30Most differential co-expression analyses rely on differential clustering; they identify clusters that contain different genes or behave differently under changing conditions or phenotypes. The most frequently used programs for differential clustering analysis, which have also been compared with others programs, are WGCNA, DICER and DiffCoEx, all of which first identify modules co-expressed across the full set of study samples. These co-expressed modules can then be correlated to predefined sample subpopulations representing, for example, disease status or tissue type.

31WGCNA determines the activity and importance of each module in each subpopulation of samples (Figure 3A and 3C). For each module, an eigengene is calculated, which is the vector that best describes the expression behaviour (in a linear fashion) of all genes within this module in the samples included in the analysis. It then prioritizes which genes in these modules are likely to underlie the phenotype associated with the module by identifying either genes behaving similarly to the eigengene of the module or those genes that are intra-modular hub genes (these tend to coincide). By design, DICER is tailored to identify module pairs that correlate differently between sample groups, e.g. modules that form one large interconnected module in one group compared with several smaller modules in another (Figure 3D). DICER may be particularly useful for time series experiments in which co-expression changes are gradual, e.g. cell cycle series experiments, where modules are specific to a particular phase and co-expressed in transitions between phases. DiffCoEx focuses on modules that are differentially co-expressed with the same sets of genes. The most extreme case of this behaviour is sets of genes that ‘hop’ from one set of correlated genes to another in a coordinated manner (Figure 3E). In this case, DiffCoEx would cluster ‘hopping’ genes in a similar manner. DINGO is a more recent tool that works similarly to DiffCoEx by grouping genes based on how differently they behave in a particular subset of samples (representing e.g. a particular condition) from the baseline co-expression determined from all samples. These are the most likely genes to explain different phenotypes that are associated with the two different networks. Each of the methods detects specific module changes by design, but they can also detect modular changes that they were not specifically designed for and may outperform other tools in the identification of these changes.

32A number of studies have used differential co-expression network analyses to identify networks unique to specific tissues or disease states. The rapid increase in publicly available RNA-seq data and projects such as GTEx and ENCODE, which generate large-scale RNA-seq profiles, has enabled co-expression analysis within and across different tissues. The GTEx project collects and provides expression data from multiple human tissues for the study of gene expression, regulation and their relationship to genetic variation. In a study comparing RNA-seq data from 35 tissues from the GTEx data set, a tissue hierarchy was constructed based on the average gene expression in each tissue. Related tissues, such as those from different brain regions, clustered together. This hierarchy was used to construct a single combined co-expression network derived from the tissue-specific co-expression networks—a meta-network. It was then shown that in tissue-specific networks, TFs with functions specific to that tissue tend to be highly expressed together with tissue-specific genes. These genes tend to form a stronger connection with each other than with other genes, but remain at the periphery of the network (thus having low centrality), while the tissue-specific TFs become more central to that module. Thus, tissue-specific TFs could be uncovered by identifying modules with increased co-expression strength in tissue-specific networks (Figure 3A and 3C) and by pinpointing the central hubs of these modules. In contrast, genes that are not TFs but are tissue-specific should be detectable by identifying genes that are at the periphery in these modules (Figure 3B). Moreover, some TFs have different roles in different tissues. These TFs would be expected to be hub genes that are central to one module under one condition and central to another module in another condition.

33Differentially connected genes are those with different co-expression partners between two sample groups. These genes appear to play a regulatory part in the difference in the phenotype observed between two groups (Figure 3D). For example, one study compared co-expression in mutant cattle with increased muscle growth with co-expression in non-mutants, using a method similar to DiffCoEx. By identifying the most differentially expressed genes and TFs showing the highest differential connection to these genes (Figure 3D), the TF containing the causal mutation (myostatin) was identified. Interestingly, the Mstn gene, which encodes this TF, hardly changed in expression itself, providing an example of how differential co-expression analysis can uncover biologically important findings not revealed by differential expression analysis alone.

34Not all methods construct a co-expression network to assess differential expression. GSNCA can be used to identify differentially co-expressed gene sets, which have to be defined a priori, between two sample groups. In the first step this method determines weight vectors for each sample group, from a correlation network. These weight vectors represent the cross-correlation of each gene with all the other genes, effectively summarizing a correlation matrix into a single vector, describing a weight for each gene. These weights for the genes representing a certain gene set are then compared between two sample groups, to determine whether the gene set is differentially co-expressed.

35Generalized Single Value Decomposition (GSVD) is a unique type of differential co-expression analysis that relies on spectral decomposition to identify modules of co-regulated genes. Unique to this approach is that it summarizes the expression of samples and all genes into a smaller number of variables, aiming to explain as much expression variation in as few variables as possible. Here we focus on the summary of gene expression into principal components or ‘genelets’, a term introduced in that can be interpreted as an analogy to co-expressed modules, and which represent the partial expression of multiple genes. The relative significance of these genelets—describing the extent to which a signal from the genelet is present (that is, the extent to which the genelet is expressed) in a data set—can be compared between two data sets. If the significance is similar, the genelet represents a co-expression pattern shared between the two data sets, whereas differences in significance indicate that the co-expression pattern is unique to one of the data sets. Higher Order (HO)-GSVD was more recently developed and uses a similar approach for comparisons between more than two data matrices.

36GSVD was first used in 2003 to analyse microarray expression data from human and budding yeast to identify common and unique pheromone and stress response patterns between these two species. HO-GSVD recently proved effective at identifying pathways important for self-renewal of neural progenitors. GSVD was shown to identify patterns unique to glioblastoma multiforme, a type of brain tumour, which was useful for prognostic purposes. Similarly, genelets that are active in normal samples were identified. These genelet signals were then removed from the total signal in cancer samples, revealing a cancer-specific signature. Both of these studies demonstrated that signatures unique to the cancer had a strong signal for genes duplicated in the cancer, as is common in cancers, suggesting that identified profiles reflect the oncogenic events in the genome.

37It is not surprising that differential co-expression methods are growing in popularity as the cost of high-quality expression data decreases. While these methods have not yet been applied to RNA-seq data, recent findings from microarray studies make this an exciting prospect. However, because these methods are sensitive to outliers, they require high-quality data.

38An alternative method for detecting differentially expressed clusters between subpopulations of data is biclustering. If a data set contains several biologically distinct but unknown sample groups, biclustering can identify genes with a similar expression pattern in only a sub-set of the samples without the need for prior sample classification (Figure 3F). This is particularly useful when such information is not available, as can be the case for large-scale single-cell RNA-seq experiments like those using the Drop-seq system or inDrop.

39In a clinical study it is often possible to predefine groups of healthy and diseased samples. However, the same disease can manifest through different mechanisms. This is a scenario common in cancer, where different mutations can lead to different alterations in co-expression patterns but a similar phenotype. Biclustering allows researchers to disentangle the mechanisms in the cases where predefining biologically relevant sample groups is difficult. For this purpose, biclustering is more effective than other co-expression analysis methods.

40Cheng et al. were first to use biclustering in co-expression analysis, followed by the development and application of many more biclustering approaches (reviewed by Pontes et al.). The choice of biclustering method depends on the number of samples and factors such as whether the samples are species- or tissue-specific and whether the included samples constitute disease phenotypes and/or different time points. Biclustering methods can be computationally challenging depending on the method used. Methods should be selected carefully because different biclustering approaches can have varying results in the same data set.

41Biclustering approaches were recently applied to RNA-seq-based expression data. Analysis of the expression data from several developmental stages of worm and fruit fly, by identifying biclusters containing similar orthologous gene sets unique to different developmental stages between the two species, led to the identification of genes with a similar, and thus conserved, function in development. Biclustering has also been applied to single-cell RNA-seq data. Because biclustering groups genes and samples simultaneously, it enabled the simultaneous identification of groups of cell types and corresponding gene modules to reveal 49 different cell types and their corresponding cell-type-specific gene modules, results that were later supported by experimental validation. With the emergence of single-cell RNA-seq, biclustering methods may be able to identify cell-type-specific modules that are present in diseased but not in healthy cells.

42Another biclustering method identified miRNAs deregulated in breast cancer through their presence in biclusters unique to cancer samples. These miRNAs have been suggested as markers for diagnosis and treatment response. Biclustering has also been used to identify tightly co-expressed sets of protein-coding genes unique to subpopulations of cancer patients, which could be used to understand patient prognosis and to further precision medicine approaches. In another cancer data set, a three-dimensional clustering method (triclustering) was used to identify genes co-expressed across subpopulations of samples and time points. This method effectively identified several known breast cancer genes in a breast cancer cell line by identifying hub genes in triclusters differentially expressed between cancer samples at early and late time points using the eigengene changes between the samples of each tricluster.

43While a comprehensive and unbiased comparison of methods used in differential co-expression analysis is desirable, the performance of the tools may be situation-dependent, varying between species, disease states and perhaps even data sets, thus making it difficult to identify the optimal method in each circumstance. An attempt was recently made to compare 10 differential co-expression algorithms, but concluded that it remains difficult to evaluate these owing to the lack of gold standard gene sets to validate the outcome of these methods. Several of the tools described in this review have been compared in publications introducing a competing method. DICER has been argued to perform better than DiffCoEx and CoXpress based on functional enrichment analysis of differentially expressed modules. HO-GSVD outperformed WGCNA and DiffCoEx based on its ability to detect clusters in simulated data. Although biclustering is a powerful approach, it does not necessarily perform better than other network analysis methods such as WGCNA, as shown by a comparison using different tools on simulated data. However, as discussed earlier, biclustering can be performed without the need for prior sample group classification.

44Although many of the tools and methods described in this review were originally created for microarray data, they are also applicable to RNA-seq data. There are RNA-seq-specific differential co-expression analysis methods, harnessing the exon- or isoform-specific expression information or allele-specific expression effects, that have been reported to perform better than other tools where this information is not considered. However, it is unclear whether these differential co-expression methods also perform better if other methods are supplied with the same isoform-specific expression information, which could be determined before differential co-expression analysis. As a result, it remains difficult to assess whether these new tools perform better than already well-established tools such as WGCNA, which can also be used on isoform-specific expression data.

45Because the tools described in this review create modules based on different criteria, it is also questionable whether the measures used in the comparisons represent desirable properties for all cases. For example, DiffCoEx groups genes based on their differential co-expression behaviour, whereas WGCNA identifies modules that are co-expressed in multiple samples and conditions. In a homogenous data set, DiffCoEx will likely detect fewer and smaller modules, indicating that there are not many differentially co-expressed genes. This will likely lead to lower enrichment scores when the performance of DiffCoEx is compared with WGCNA on such homogenous data, whereas it is merely an indication that not many co-expression partner changes occur in the data. By contrast, tools that focus on shared co-expression modules are likely to find strong correlation modules with high enrichment scores, which may not be relevant if the goal is to identify regulatory modules. WGCNA has been widely shown to perform well under many different circumstances and for different purposes. However, it requires information on the sample conditions to assign modules to conditions. If this information is not available (as in large-scale single-cell RNA-seq experiments) or if researchers wish to identify subgroups within the sample groups, biclustering is a more suitable approach.

46To get a systematic assessment of the performance of different tools and methods, projects such as DREAM4 and DREAM5 have been invaluable. These challenged researchers to construct regulatory networks from simulated and in vivo benchmark data sets. As these challenges are predefined they allow researchers to test their methods/tools in an unbiased manner. However, these challenges were last posed in 2010 and many new methods and tools have been developed since.

47Experimental validation often focuses on single genes. As these experiments are costly and time-consuming, high confidence predictions of causal genes are of great importance. An analysis based solely on co-expression does not (yet) provide this level of confidence. Therefore, incorporation of information from other types of data can help to prioritize genes that may underlie a phenotype. This can be achieved, for example, using information describing which genes are TFs, as is done for regulatory predictions by GENIE3. However, a focus on TFs is rarely sufficient, and integration of multiple data types is often required to increase the accuracy and usefulness of the resulting networks.

48Genome-wide transcription factor binding site (TFBS) analysis was introduced in the beginning of this millennium using chromatin immunoprecipitation followed by microarray analysis, also known as ChIP-chip, which was later replaced by the more accurate ChIP-seq. These data were used to create a genome-wide integrated regulatory network from gene expression and TFBS data. Combined analysis of ChIP-chip-based TFBSs and expression data initially showed that, in 58% of the cases, the TFs bound to the promoter region of the gene were indeed regulated by the corresponding TF. A partial least squares approach (a well-known method for analysis of high-dimensional data with several continuous response variables) was later proposed to identify false positives and distinguish the activation and repression activities of TFs. A more recent method harnesses the rapidly increasing availability of ChIP-seq data in combination with expression data to rank the genes bound by a TF, which can be used to prioritize the most likely TF targets. Tools to conduct similar analyses, integrating expression and ChIP data, have also been published.

49Independent from the approach used to identify them, network modules can be further investigated for shared eQTL gene targets, TF/miRNA targets or enriched binding motifs. Several computational methods and publicly available data sets are available for multi-omics data integration. For example, information about eQTLs can be acquired from recent large-scale blood-based trans-eQTL meta-analysis or eQTL studies conducted in other tissue types. Transcription factor binding sites (TFBSs) can be collected from databases such as JASPAR and DeepBind, which consist of TF binding motifs inferred from experimental data. Binding sites can be further prioritized by investigating tissue-specific ChIP-seq peaks from ENCODE. Finally, miRNA–target interactions can be identified using several in silico target prediction tools or using manually curated databases of experimentally supported target interactions.

50Combining information from different layers of data may lead to new biologically interpretable associations in a number of ways. If intra-modular hub genes are TFs or targets of a TF, this TF is more likely to have a causal role in the phenotype under investigation. If multiple Genome-Wide Association Study (GWAS) hits exist in the same module, their cumulative presence can significantly contribute to disease development. Differential methylation states of genes within a co-expression module can elucidate methylation patterns underlying disease. If multiple genes are regulated by the same genetic variant (under a trans-eQTL effect), it may be possible to identify the gene responsible for the alterations of the network by identifying the cis-eQTL gene driving the trans-eQTL effects (Figure 4). This is supported by the fact that genes under trans-regulation of disease-associated genomic variants are sometimes functionally connected with the processes or pathways associated with the corresponding disease. Good examples of this are IFN (interferon)-α and complement pathways in which several genes were under trans-regulation of a systemic lupus erythematosus-associated variant, possibly via cis-regulation of IKZF1. The integration of regulatory genetic variant information into co-expression network analysis, with cis-eQTLs used as causal anchors, identified TYROBP as the most likely causal factor in late-onset Alzheimer disease patients, a finding supported by the observation that mutations in this gene are known to cause Nasu-Hakola disease. Lastly, copy number variation can affect gene expression levels, and including such information may help identify and/or explain alterations in co-expression network structures present in diseases or traits.

51Overall, integration of multiple data types increases the accuracy of the resulting predictions. For example, modules unique to different subtypes of cancer were identified by integrating tumour genome sequences with gene networks, and these modules may be useful for prognosis and identification of putative targets for personalized medicine-based treatments. A number of tools, described earlier in this review, can be used for differential co-expression analysis, but can also be applied to other data types. In the initial DINGO publication, the authors conducted a combined analysis on mRNA expression, DNA copy number variation and methylation data. By overlaying the differential networks of each data type and identifying edges present in all of them, a number of genes from the PI3K pathway were identified as important players in glioblastoma multiforme patients. This pathway is an already-established therapeutic target, supporting the notion that this is an effective approach for identifying relevant targets for disease studies. A recently published tool, CoRegNet, allows the integration of different types of data in a co-expression analysis by identifying co-operative regulators of genes from different data types. Another established approach, cMonkey, achieves similar data integration by calculating the joint bicluster membership probability from different data types by identifying groups of genes that group together in multiple data types.

52In recent years, differential co-expression analyses have been increasingly used to analyse large data sets. This may be attributed to the decreased costs of large-scale gene expression profiling, in particular RNA-seq, to increased sample sizes, and to the greater availability of tissue-specific data from perturbation experiments, which are required for fruitful differential co-expression analyses. Likewise, biclustering algorithms have benefitted from larger sample sizes and higher data quality, as shown by the identification of co-expressed modules unique to cancer subtypes. The usefulness of biclustering on single-cell RNA-seq data has been demonstrated by the classification of different cell types and by the identification of clusters of genes uniquely co-expressed in specific cell types. We expect these approaches to be more widely applied in the future, as they benefit from an increase in RNA-seq data quantity and quality, which will allow for more accurate identification of tissue-specific and cell-type-specific disease-related modules and regulators.

53Large-scale single-cell sequencing technology is increasingly used and the first co-expression studies using such techniques have uncovered cell-type-specific co-expression modules that would have gone undetected in multi-cell-type co-expression analyses. Because the latter represent the aggregated signals of multiple cell types, they usually cannot detect alterations in cell subpopulations between different experimental groups. This is supported by the observation that the expression of cell cycle genes associated with ageing decreased in the analysis of non-cell-type-specific data. However, data from single-cell experiments revealed that this observation was caused by a decreased proportion of the G1/S cells that highly express cell cycle genes rather than by altered expression across the whole cell population.

54An additional prospect is the detection of mutations from RNA-seq data. As mutations accumulate with age in different cells, these can be used to identify the origin of the cell. Mutation accumulation has been used to study cancer development and the origin of metastases. In large-scale single-cell RNA-seq experiments, mutations could be used to separate cells based on their origin, or to group cells based on the mutations they harbour. Cells harbouring the same mutations can be investigated for co-expression patterns, and modules unique to cells with a specific mutation may be detected. This may allow the direct linking of mutations to expression modules, with the limitation that only mutations in coding regions are detectable in RNA-seq data.

55Although there are many exciting new possibilities with single-cell RNA-seq data, important challenges remain. Typically, a low number of reads per cell are sequenced and then the signal from multiple cells of the same type is aggregated to acquire a cell-type-specific gene-expression profile. It is hard to acquire sufficient data for rarer cell populations, such as stem cells, and this is currently limiting analyses on these cell types. Additionally, the low number of reads per cell leads to sparse expression matrixes to which normalization methods currently used in canonical RNA-seq analyses are not attuned. These normalization methods often also assume that the majority of genes do not change in expression between different samples, which is not necessarily the case in single-cell RNA-seq owing to variation in expression across different cells. This is further exacerbated by the difficulty in obtaining high-quality RNA from single cells. These and other issues are further discussed in.

56In addition to the normalization issues that occur in single-cell RNA-seq, the optimal method for normalizing bulk RNA-seq data is also still not clear. The widely used Fragment/Reads Per Kilobase Million (FPKM) normalization has been debated and although alternatives have and are being created, each method has its limitations. Additionally, from our experience, the use of different mapping tools can in some cases lead to different results. Although some comparisons between different tools and methods have been made, a large-scale comparison, using e.g. public data, would identify such cases and define best practices for pursuing each research question.

57With the increased availability of different data types such as RNA-seq, genome sequences, ChIP-seq, methylome and proteome data, it will become possible to integrate these data sets to more accurately predict regulatory genes. Projects from large consortia like GTEx, the Epigenome Roadmap and ENCODE are already generating data from multiple-omics levels that facilitate these integrated analyses. To identify regulatory relationships, perturbation data are preferable, as canonical data cannot distinguish between true and false positives in regulatory relationships. Furthermore, regulatory relationships can be highly cell-type-, tissue- or developmental-stage-specific. Only a handful of tools and methods are currently available to investigate multi-omics data, and the tools that exist mostly integrate only two layers of omics data. Integrated network analyses come with additional mathematical challenges, and best practices are far from established. Further research on this topic is of great interest to the research community, as it will allow a better understanding of regulatory mechanisms that can explain co-expression patterns and disease mechanisms. A better understanding of these disease mechanisms and corresponding co-expression patterns will facilitate the identification of appropriate targets for intervention studies.

58RNA-seq-based co-expression analysis can be used to assign putative functions to non-coding RNAs and to identify candidates for roles in disease.In co-expression networks, hub gene identification has a limited power for identifying targets for follow-up studies; yet, this can be enhanced by integrated network analyses, which may incorporate GWAS hits, eQTLs, TFBSs and other data layers.Differential co-expression analyses can reveal genes that have different co-expression partners between healthy and disease state and can help to uncover regulators underlying disease and other phenotypes.Methods such as biclustering and Generalised Single Value Decomposition (GSVD) allow the identification of signals/modules unique to specific cancer subtypes, which may serve a purpose in prognosis and for precision medicine.

59RNA-seq-based co-expression analysis can be used to assign putative functions to non-coding RNAs and to identify candidates for roles in disease.

60In co-expression networks, hub gene identification has a limited power for identifying targets for follow-up studies; yet, this can be enhanced by integrated network analyses, which may incorporate GWAS hits, eQTLs, TFBSs and other data layers.

61Differential co-expression analyses can reveal genes that have different co-expression partners between healthy and disease state and can help to uncover regulators underlying disease and other phenotypes.

62Methods such as biclustering and Generalised Single Value Decomposition (GSVD) allow the identification of signals/modules unique to specific cancer subtypes, which may serve a purpose in prognosis and for precision medicine.