--- title: "Gene Co-Expression Analysis" canonical: "https://help.biobam.com/space/OED0324/4935614466/Gene+Co-Expression+Analysis" format: markdown --- ## Introduction This tutorial provides a step-by-step guide to performing a [gene co-expression analysis with WGCNA](https://docs.omicsbox.biobam.com/latest/Gene-Coexpression-WGCNA/) within OmicsBox. WGCNA (Weighted Gene Co-expression Network Analysis) groups genes into modules of highly co-expressed genes, relates these modules to sample traits, and identifies the most connected genes (hub genes) within each module. This makes it well suited to reduce a large expression matrix to a set of biologically meaningful modules and to identify gene groups and candidate regulators associated with a condition of interest. Please cite WGCNA as: > Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008 Dec 29;9:559. doi:10.1186/1471-2105-9-559 > > Zhang B, Horvath S. A general framework for weighted gene co-expression network analysis. Stat Appl Genet Mol Biol. 2005;4:Article17. doi:10.2202/1544-6115.1128 ### **Dataset** The dataset used in this tutorial is derived from a microarray study of the stepwise progression of HCV-related hepatocellular carcinoma (HCC). The design is as follows: - **Sample type**: human liver tissue. - **Experimental conditions**: eight consecutive stages of hepatocarcinogenesis, from pre-neoplastic lesions to end-stage tumor: normal, cirrhosis without HCC, cirrhosis, low-grade dysplastic, high-grade dysplastic, very early and early HCC, advanced HCC, and very advanced HCC. - **Total Samples**: 75 tissue samples across the eight stages. - **Platform**: Affymetrix Human Genome U133 Plus 2.0 array. In this tutorial, the pathological stage is used as the trait against which co-expression modules are correlated, allowing modules and their hub genes to be linked to specific points in disease progression. > Data Source: The expression data were originally published by [Wurmbach E, et al. (2007)](https://doi.org/10.1002/hep.21622) in "Genome-wide molecular profiles of HCV-induced dysplasia and hepatocellular carcinoma" (Hepatology 45(4):938–947), and are publicly available in the NCBI GEO database under accession number [GSE6764](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE6764). > > The WGCNA workflow reproduced here follows the analysis of [Yin L, et al. (2018)](https://doi.org/10.3390/genes9020092), "Identification of Key Pathways and Genes in the Dynamic Progression of HCC Based on WGCNA" (Genes 9(2):92). ## Bioinformatic Analysis The analysis pipeline starts from raw sequencing reads, which first need to be cleaned and then quantified to obtain an estimate of the gene or transcript expression for each sample. These first steps are covered in the [Gene-level](https://docs.omicsbox.biobam.com/latest/Gene-level-Quantification/) and [Transcript-level](https://docs.omicsbox.biobam.com/latest/Transcript-level-Quantification/) Analysis Example Datasets, so this tutorial will start directly from a count table. Microarray data can also be used as input, as is the case for the dataset in this tutorial; in that case the expression matrix is not generated within OmicsBox and must be [loaded](https://docs.omicsbox.biobam.com/latest/Load/#load-rna-seq-count-table) directly. Then, WGCNA groups genes into modules of highly co-expressed genes and relates these modules to the sample traits, here the pathological stages of HCC progression, highlighting the modules and hub genes most associated with the disease. Finally, the [functional enrichment analysis](https://docs.omicsbox.biobam.com/latest/Gene-Coexpression-WGCNA/#fishers-exact-test) provides biological meaning to those modules by identifying the biological processes and pathways they represent. ### **1. Renaming Input Data** The count table for this dataset uses Affymetrix probe set IDs as feature names, which are not intuitive to work with or interpret. In this first step, these probe set IDs are mapped to gene accession IDs using the platform annotation, so that modules and hub genes can be interpreted directly in the following sections. First, load the count table. Go to Transcriptomics > Load > [Load RNA-Seq Count Table](https://docs.omicsbox.biobam.com/latest/Load/#load-rna-seq-count-table) and select the input file in txt format. Leave all parameters by default. Once loaded, in the count table Side Panel, go to Actions > Rename Features. #### Input - [Count Table](https://drive.google.com/file/d/18GQT9Rc57tFXZORWOzcAWTr834hRZ6sY/view?usp=sharing): microarray gene expression values across all 75 samples. - [Mapping File](https://drive.google.com/file/d/1LMg-kjQ-wgBVb6oP8JAfjYkEmd93vHfI/view?usp=sharing): a two-column text file mapping each Affymetrix probe set ID to its corresponding gene accession ID. It is derived from the platform annotation for the Affymetrix Human Genome U133 Plus 2.0 array ([GPL570](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL570)), keeping only the first two columns of the file downloaded from NCBI. #### Parameters - Option: Mapping - Aggregation Method: Sum Rows - Mapping [File: provid](#)e here the input mapping file. All other parameters are left at their default values. #### Output - [Renamed Count Table](https://drive.google.com/file/d/1e1THNIkK2srIoFRRjmwHdAb8x_Wu7vTU/view?usp=sharing): a new count table, opened in a separate tab, in which features are labeled with gene accession IDs. Because several probe sets can map to the same gene, their values are combined using the Sum Rows aggregation method. This renamed table is the input for the co-expression analysis. ![image-20260706-135354.png](media://0c2a5c95-485e-4d68-8887-ef2f47725950) ### 2. Gene Co-Expression Analysis To run the analysis, open the renamed count table and go to the Side Panel > Actions > [Gene Co-Expression Analysis](https://docs.omicsbox.biobam.com/latest/Gene-Coexpression-WGCNA/). #### Input - [Renamed Count Table](https://drive.google.com/file/d/1e1THNIkK2srIoFRRjmwHdAb8x_Wu7vTU/view?usp=sharing): renamed microarray gene expression values across all 75 samples. - [Experimental Design](https://drive.google.com/file/d/1oRJmHO46gBmqk-eVNJtHW3X5S9IzqYov/view?usp=sharing): structured text file mapping each sample to its pathological stage, used later to correlate the modules with disease progression. #### Parameters - High Variable Genes Method: Top N Genes - N. High Variable Genes: 17000 - Power Calculation: Manual - Power for Soft Threshold: 8 - Merge Cut Height: 0.2 All other parameters are left at their default values. These settings follow those used in the original study by Yin et al. (2018). #### Output - [Gene Co-Expression Results](https://drive.google.com/file/d/1LgexvVNP2XfVOeYASY1EEgp2v5qGkJT2/view?usp=sharing): this is the primary output of the analysis, which contains the statistics per gene. Each gene is assigned to a module (identified by a color), and, optionally, to a hub tag. The hub tag marks the genes with the highest connectivity within their module, that is, those whose expression profile correlates most strongly with the module eigengene; these hub genes are often the most biologically relevant candidates. When an experimental design is provided, as in this case, the table also includes, for each gene, the correlation and its associated p-value with each experimental condition. The table can be filtered and sorted by any column by clicking on its header. The Side Panel shows the available charts and downstream analyses. - [Summary Report](https://drive.google.com/file/d/1y1RpiBytR9b5L3Oj3fY1X9fhKqQLbogh/view?usp=sharing): it contains basic statistics about the run, including the total number of identified modules and the number of filtered genes and samples, if any. It also includes general statistics per module, where the module sizes and the average module membership (Avg MM) can be seen. In this analysis, 22 modules were generated, ranging in size from 42 to 3,065 genes, with Avg MM values between 0.6 and 0.8, indicating mid-defined, internally coherent modules. ![image-20260706-134755.png](media://76e7d89e-234a-4f4c-b45a-ca0509839425) ##### Assessing the Parameter Choices **Gene Variance Elbow Plot** This plot ranks genes by decreasing variability, with genes ordered from highest to lowest variability on the x-axis and their variability values on the y-axis. Retained genes are shown in green and excluded genes in grey. The curve usually declines steeply and then flattens; the "elbow" marks the transition between the two regions. Genes before the elbow carry meaningful biological variation, whereas genes beyond it show similar, low variability and add little to the analysis. In this case, the 17,000 genes retained fall on the flat part of the curve, indicating that many of them contribute little additional information. A value closer to the elbow, around 10,000 genes, would provide a more appropriate cut-off while retaining the most informative genes. ![image-20260706-102000.png](media://5af66425-ff25-45f6-8876-8a03c2f71419) **Soft Threshold Power plot** This plot supports the selection of the soft-thresholding power used for network construction. For each tested power, it reports the Scale-Free Topology Fit Index (R², left y-axis) and the Mean Connectivity (right y-axis). A higher R² indicates a better fit to the scale-free topology expected in biological networks, where a few hub genes show many connections and most genes show few. As the power increases, weak correlations are suppressed and the Mean Connectivity decreases. A suitable power balances a high R² against a reasonable connectivity, and is usually the lowest power at which R² reaches a high level (commonly ≥ 0.85, marked by the dashed grey line) or begins to plateau. The view can be adjusted using the plot Side Panel. For instance, we can increase the “Label Interval (bottom)” parameter to see all the Power values on the x-axis. In this case, the selected power of 8 lies well above the R² threshold, at a point where the Mean Connectivity is already low. A more conservative value, around 5, would sit closer to the point where the scale-free topology criterion is first met while preserving greater connectivity. ![image-20260706-102227.png](media://a13062ed-e6b1-4e71-9a7f-6dc789f3287e) > ℹ️ Taken together, these two diagnostic plots, along with the average module membership values not being particularly high (0.8 at most), suggest that the analysis could be re-run with adjusted parameters to obtain more clearly defined modules. For the sake of brevity, this tutorial retains the parameters used in the original study by Yin et al. (2018), but these diagnostics illustrate how the plots can guide parameter selection in other analyses. ##### Exploring the Modules Two plots help identify the most relevant modules and the relationships between them: the Module–Trait Correlation plot and the Modules Correlation plot. **Module–Trait Correlation** This plot visualizes the relationship between each module and the experimental conditions, with the traits on the y-axis and the modules on the x-axis. Each bubble represents the correlation between a module and a trait: the color indicates the correlation coefficient (red for positive, blue for negative, near-white for weak or no correlation), and the bubble size reflects its magnitude, computed from the p-value through a Sigmoid function. It is therefore a useful starting point to identify modules worth studying further. The minimum and maximum bubble sizes have been set to 5 and 30, respectively, though the plot Side Panel, which produces a cleanier view of the bubbles. Other values may be more appropiate depending on the display. Here, several modules are associated with very advanced HCC, most notably the turquoise module. Some of these are also associated with early and advanced HCC, whereas the yellow and royalblue modules are associated with very advanced HCC only. Additional sets of modules are associated with the low-grade dysplastic and cirrhotic stages. These stage-specific associations point to modules whose co-expressed genes may be involved in particular steps of disease progression. ![image-20260706-113312.png](media://5737b21e-8f16-4973-8e57-292c2cdf03ef) **Modules Correlation** This plot shows the direction and strength of the correlations between modules, with the identified modules on both axes. As above, the color indicates the correlation coefficient and the bubble size reflects the statistical significance, computed from the p-value through a Sigmoid function. Only one half of the matrix is shown, with redundant pairs omitted. Positive correlations are generally the most straightforward to interpret, as positively correlated modules tend to share expression patterns and may reflect coordinated biological activity or common regulation. The minimum and maximum bubble sizes have been set to 5 and 25, respectively, trought the plot Side Panel, which produces a cleanier view of the bubbles. Other values may be more appropiate depending on the display. These results are consistent with the module–trait associations: modules associated with the same trait tend to be correlated with each other. The turquoise module, for example, is correlated with the other modules linked to advanced disease, though additional correlations between modules can also be observed. ![image-20260706-113605.png](media://c2964375-b34e-400a-b191-df92ce15dc57) ##### Visualizing Networks **Module(s) Network** The Module(s) Network plot displays the co-expression network of the genes belonging to one or more selected modules. Hub genes are shown as diamonds and non-hub genes as circles, node size reflects the Module Membership (MM), edge width represents the connection strength between genes, and node color indicates the module. The genes included can be controlled with the Module(s) to Plot, Min. Connectivity, and Min. MM parameters, and the appearance can be further adjusted from the network Side Panel. In this case, the turquoise module was selected for visualization. As its average module membership was somewhat low, the Min. MM threshold was lowered to 0.8 so that the principal genes of the module are displayed. The resulting network shows these genes and the strength of their connections through the edges. If a particular gene of interest is identified, it can also be examined individually using the Gene Network plot. These visualizations are intended as an exploratory aid to inspect module structure and identify well-connected genes. ![image-20260706-115827.png](media://de2b30a0-147a-4023-8574-3817450bbcd3) ### 3. Downloading Functional Annotation from BioMart To relate the module genes to biological functions, a reference associating each gene with its function is needed. This can be generated from scratch using the tools in the [Functional Analysis module](https://docs.omicsbox.biobam.com/latest/Module-Functional-Analysis/) (Blast2GO, InterProScan, eggNOG), or an existing annotation can be retrieved from BioMart, as done here. Go to Functional Analysis > Load > [Load data from BioMart](https://docs.omicsbox.biobam.com/latest/Load-Data/#load-data-from-biomart). #### Parameters - Origin: Vertebrates - Dataset (Mart): Ensembl Genes 116 - Database: Human genes (GRCh38.p14) - Type of ID: AFFY HG U133 Plus 2 probe ID(s) [e.g. 1553551_s_at] - Retrieve Sequence Data: No - Which IDs to Import: Whole Genome All other parameters are left at their default values. #### Output - [Functional Annotation](https://drive.google.com/file/d/1o_7aiDqHMUGyVbuKfqI53zmQW1rYfiay/view?usp=sharing). In this project, the Affymetrix probe set IDs are associated with their Gene Ontology (GO) terms and Enzyme Codes (EC). As with the count table, the features are labeled with Affymetrix probe set IDs and must be **renamed** to match the rest of the data. In the Side Panel, go to Tools > [Batch Rename](https://docs.omicsbox.biobam.com/latest/Tools/#batch-rename) and provide the same mapping file used previously. #### Parameters - Option: Replace - Use Mapping [File: checke](#)d - Mapping File: [select the mapping file](https://drive.google.com/file/d/1LMg-kjQ-wgBVb6oP8JAfjYkEmd93vHfI/view?usp=sharing). All other parameters are left at their default values. #### Output - [Renamed Functional Annotation](https://drive.google.com/file/d/1o_7aiDqHMUGyVbuKfqI53zmQW1rYfiay/view?usp=sharing): a project associating each gene, now labeled with its gene accession ID, with its GO terms and Enzyme Codes. This annotation is the functional reference used to interpret the modules, so it should be saved and kept available for the enrichment analysis in the following step. ![image-20260706-140337.png](media://e6cdc4ef-5dff-4fa6-9e8d-1f11b314ff90) ### 4. Functional Enrichment Functional enrichment analysis identifies which biological functions are statistically overrepresented in a specific group of genes compared to a reference background. This step helps translate the genes grouped in a module into meaningful biological processes, molecular functions, or cellular components. In OmicsBox, this analysis is performed using a Fisher's Exact Test, which can be launched from the Gene Co-Expression results Side Panel > Actions > [Fisher's Exact Test](https://docs.omicsbox.biobam.com/latest/Gene-Coexpression-WGCNA/#fishers-exact-test). #### Input - [Functional Annotation](#), generated in step 3. A file containing all the genes along with their annotated biological functions, such as Gene Ontology (GO) and Enzyme Commission (EC) terms. #### Parameters - Test Module(s): yellow, royalblue - Test Hub Features Only: checked - Reference Annotation: select the file provided in the input section. All other parameters are left at their default values. Here, the yellow and royalblue modules are tested, as both were associated specifically with the very advanced HCC stage. #### Output - [Fisher's Exact Test Results](https://drive.google.com/file/d/1hsyl80Pk0ifNOU6lrwhi3NGDBKhpAonW/view?usp=sharing). Table listing all the functional terms tested, with a Tag specifying whether each term is statistically over- or under-represented in the analyzed modules. To obtain more concrete results, the enriched terms can be summarized. From the Fisher's Exact Test results, go to the Side Panel > Actions > Reduce to Most Specific. This removes redundant parent terms and keeps the most specific ones, making the results easier to interpret. ![image-20260706-134358.png](media://aac31dd3-0c22-4f3a-be20-050ecbeafe64) The same procedure can be applied to any other module of interest, for instance those associated with the low-grade dysplastic or cirrhotic stages, to explore the biological functions characterizing each step of disease progression. ## **CloudUnits Summary** OmicsBox utilizes [OmicsCloud](https://www.biobam.com/omicscloud-cloud-platform/), an AWS-based high-performance computing (HPC) infrastructure, to execute resource-intensive bioinformatics tasks. These costs are measured in [Cloud Units](https://www.biobam.com/cloud-computation/), which directly correspond to the compute infrastructure used, such as CPU seconds and data transfer volume. Most features within OmicsBox are included for free with an active subscription; only features related to Sequence Alignment and Assembly consume Cloud Units. Different tools need more or less computational resources, depending on the algorithm, input data, and used parameters. The following table summarizes the resources consumed during this study. Please note that these are approximate values intended to serve as a reference: | | | | | | | --- | --- | --- | --- | --- | | **Tool** | **Input Size** | **Execution Time** | **Cloud Units Consumed** | **Cloud Units Charged** | | Gene Co-Expression Analysis | 11.2 MiB | 26m | 2424 | 0 | The BioMart data download and Fisher’s Exact Test do not run in the OmicsCloud.