Application of weighted gene co-expression network analysis to identify key modules and hub genes in oral squamous cell carcinoma tumorigenesis

Figures & data

Figure 1 Clustering of samples and determination of soft-thresholding power.

Notes: (A) The clustering was based on the expression data of GSE30784, which contained 167 OSCC, 17 dysplasia and 45 normal samples. The top 5,000 genes with the highest SD values were used for the analysis by WGCNA. The color intensity was proportional to disease status (normal, dysplasia and OSCC), sex and age. (B) Analysis of the scale-free fit index for various soft-thresholding powers (β). (C) Analysis of the mean connectivity for various soft-thresholding powers. In all, 9 was the most fit power value.
Abbreviations: OSCC, oral squamous cell carcinoma; WGCNA, weighted gene co-expression network analysis.

Figure 2 Construction of co-expression modules by WGCNA package in R.

Notes: (A) The cluster dendrogram of module eigengenes. (B) The cluster dendrogram of genes in GSE39784. Each branch in the figure represents one gene, and every color below represents one co-expression module.
Abbreviation: WGCNA, weighted gene co-expression network analysis.

Figure 3 (A) Interaction relationship analysis of co-expression genes. Different colors of horizontal axis and vertical axis represent different modules. The brightness of yellow in the middle represents the degree of connectivity of different modules. There was no significant difference in interactions among different modules, indicating a high-scale independence degree among these modules. (B) Heatmap of the correlation between module eigengenes and the disease status of OSCC. The turquoise module was the most positively correlated with status, and the brown module was the most negatively correlated with status. (C) Hierarchical clustering of module hub genes that summarize the modules yielded in the clustering analysis. (D) Heatmap plot of the adjacencies in the hub gene network.

Abbreviation: OSCC, oral squamous cell carcinoma.

Figure 4 (A) Scatter plot of module eigengenes in the turquoise module. (B) Scatter plot of module eigengenes in the brown module. (C) The top 30 hub genes in the turquoise module. (D) The top 30 hub genes in the brown module. Nodes represent genes, and node size is correlated with connectivity of the gene by degree.

Figure 5 Survival analysis of hub genes.

Notes: (A–E) Five hub genes in the turquoise module with the highest node degree and significant results of survival analysis (P⩽0.05 was regarded as significant). They were MMP1, TNFRSF12A, PLAU, FSCN1 and PDPN, respectively. (F–J) Five hub genes in the brown module with the highest node degree and significant results of survival analysis (P⩽0.05 was regarded as significant). They were KRT78, EVPL, GGT6, SMIM5 and CYSRT1, respectively.

Table 1 GO enrichment analysis of turquoise module and brown module

Category	Term	Count	%	P-value
Brown module
BP	GO:0055114~oxidation–reduction process	66	7.03	2.20E–10
BP	GO:0008544~epidermis development	19	2.02	7.42E–08
BP	GO:0030855~epithelial cell differentiation	14	1.49	2.17E–05
BP	GO:0098609~cell–cell adhesion	29	3.09	9.53E–05
BP	GO:0006805~xenobiotic metabolic process	13	1.38	2.92E–04
BP	GO:0018149~peptide cross-linking	10	1.06	5.07E–04
BP	GO:0030216~keratinocyte differentiation	12	1.27	8.76E–04
BP	GO:0033539~fatty acid beta-oxidation using acyl-CoA dehydrogenase	6	0.63	1.22E–03
BP	GO:0030148~sphingolipid biosynthetic process	9	0.95	1.74E–03
BP	GO:0048662~negative regulation of smooth muscle cell proliferation	7	0.74	2.09E–03
CC	GO:0070062~extracellular exosome	240	25.58	7.17E–22
CC	GO:0005789~endoplasmic reticulum membrane	78	8.31	4.60E–08
CC	GO:0016324~apical plasma membrane	36	3.83	3.62E–07
CC	GO:0030057~desmosome	9	0.95	8.84E–06
CC	GO:0005913~cell–cell adherens junction	34	3.62	2.62E–05
MF	GO:0005198~structural molecule activity	31	3.30	2.55E–06
MF	GO:0016491~oxidoreductase activity	27	2.87	3.38E–06
MF	GO:0098641~cadherin binding involved in cell–cell adhesion	31	3.30	5.91E–05
MF	GO:0005509~calcium ion binding	56	5.97	3.58E–04
MF	GO:0003824~catalytic activity	21	2.23	6.79E–04
Turquoise module
BP	GO:0050900~leukocyte migration	20	3.82	8.08E–10
BP	GO:0031581~hemidesmosome assembly	8	1.52	7.55E–09
BP	GO:0045926~negative regulation of growth	8	1.52	4.07E–07
BP	GO:0030335~positive regulation of cell migration	19	3.63	3.00E–06
BP	GO:0001525~angiogenesis	21	4.01	3.27E–06
BP	GO:0007155~cell adhesion	32	6.11	3.54E–06
BP	GO:0007160~cell–matrix adhesion	13	2.48	5.85E–06
BP	GO:0030198~extracellular matrix organization	19	3.63	7.27E–06
BP	GO:0071294~cellular response to zinc ion	7	1.33	7.99E–06
BP	GO:0001666~response to hypoxia	17	3.25	2.00E–05
CC	GO:0005925~focal adhesion	47	8.98	1.97E–17
CC	GO:0001726~ruffle	13	2.48	4.79E–06
CC	GO:0001725~stress fiber	10	1.91	1.18E–05
CC	GO:0008305~integrin complex	7	1.33	6.48E–05
CC	GO:0005938~cell cortex	13	2.48	1.14E–04
MF	GO:0005178~integrin binding	14	2.67	4.89E–06
MF	GO:0050431~TGF-beta binding	6	1.147	4.75E–05
MF	GO:0031418~L-ascorbic acid binding	6	1.14	1.98E–04
MF	GO:0001968~fibronectin binding	6	1.14	5.73E–04
MF	GO:0051015~actin filament binding	12	2.29	8.96E–04

Abbreviations: GO, Gene Ontology; BP, biological process; CC, cellular component; MF, molecular function.

Table 2 KEGG pathway enrichment analysis of turquoise module and brown module

Category	Term	Count	%	P-value
Brown module
KEGG	hsa01100: metabolic pathways	107	11.41	6.93E–10
KEGG	hsa00100: steroid biosynthesis	8	0.85	3.03E–05
KEGG	hsa00280: valine, leucine and isoleucine degradation	10	1.07	4.11E–04
KEGG	hsa01130: biosynthesis of antibiotics	22	2.35	1.78E–03
KEGG	hsa00561: glycerolipid metabolism	10	1.07	1.99E–03
KEGG	hsa03320: PPAR signaling pathway	10	1.07	5.40E–03
KEGG	hsa04071: sphingolipid signaling pathway	14	1.49	6.16E–03
KEGG	hsa00350: tyrosine metabolism	7	0.75	6.75E–03
KEGG	hsa04530: tight junction	15	1.60	7.62E–03
KEGG	hsa00600: sphingolipid metabolism	8	0.85	7.78E–03
Turquoise module
KEGG	hsa04510: focal adhesion	27	5.16	3.64E–10
KEGG	hsa04151: PI3K–Akt signaling pathway	30	5.74	3.84E–07
KEGG	hsa04512: ECM–receptor interaction	14	2.68	1.64E–06
KEGG	hsa05200: pathways in cancer	30	5.74	5.52E–06
KEGG	hsa05205: proteoglycans in cancer	20	3.82	7.96E–06
KEGG	hsa05222: small-cell lung cancer	12	2.29	4.27E–05
KEGG	hsa04978: mineral absorption	9	1.72	5.74E–05
KEGG	hsa04810: regulation of actin cytoskeleton	15	2.87	0.004426
KEGG	hsa05202: transcriptional misregulation in cancer	13	2.49	0.004612
KEGG	hsa05412: ARVC	8	1.53	0.005263

Abbreviations: KEGG, Kyoto Encyclopedia of Genes and Genomes; ECM, extracellular matrix; ARVC, arrhythmogenic right ventricular cardiomyopathy.

Figure 6 Validation of hub genes in GSE31056.

Notes: (A) Volcano plot visualizing DEGs in GSE31056 (23 normal samples and 21 OSCC). The vertical lines demarcate the fold change values. The right vertical line corresponds to more than or equal to twofold up changes and the left vertical line more than or equal to twofold down changes, while the horizontal line marks a –log₁₀ adjusted P-value of 0.01. (B) Heat map hierarchical clustering reveals DEGs in OSCC groups compared with those in control groups. (C) Identification of common genes between DEGs and the turquoise module by overlapping them. The five hub genes in the turquoise module were also DEGs in GSE31056. (D) Identification of common genes between DEGs and the brown module by overlapping them. The five hub genes in the brown module were also DEGs in GSE31056.
Abbreviations: DEGs, differently expressed genes; OSCC, oral squamous cell carcinoma.

Figure 7 Validation of hub genes in the transcriptional level.

Notes: (A) Validation of hub genes in GSE30784. The expression status of five hub genes in the turquoise module was positively correlated with the disease status, indicating they played an important role in promoting tumorigenesis. This was in accordance with the results of WGCNA. The expression status of five hub genes in the brown module was negatively correlated with the disease status, indicating they played an important role in suppressing tumorigenesis. This was in accordance with the results of WGCNA. (B) Validation of hub genes in GSE74530 and the results were the same as earlier. (C) Validation of hub genes by TCGA HNSC data in GEPIA and we obtained similar results as above. These findings indicated that our results were reliable. *** indicates ⩽0.001, **** indicates ⩽0.0001.
Abbreviations: WGCNA, weighted gene co-expression network analysis; HNSC, head and neck squamous cell carcinoma.

Figure 8 Validation of hub genes in the translational level.

Notes: (A–D) Validation of five hub genes in turquoise module by The Human Protein Atlas database (IHC). There was no related IHC samples of MMP1 in the database. The translational expression level of the rest four hub genes was positively correlated with disease status as they were upregulated in OSCC samples. (E–I) Validation of five hub genes in the brown module by The Human Protein Atlas database (IHC). The translational expression level of five hub genes was negatively correlated with disease status, as they were downregulated in OSCC samples.
Abbreviations: IHC, immunohistochemistry; OSCC, oral squamous cell carcinoma.

Figure 9 Genetic alterations associated with hub genes in TCGA HNSC.

Notes: (A) A visual summary across on a query of 10 hub genes showing genetic alteration of 10 hub genes in TCGA HNSC patients. (B) The total alteration frequency of 10 hub genes in TCGA HNSC is illustrated. (C) The network contains 55 nodes, including our 10 query genes and the 50 most frequently altered neighbor genes (only four out of 10 were correlated with the 50 genes). Relationship between hub genes and tumor drugs is also illustrated.
Abbreviation: HNSC, head and neck squamous cell carcinoma.

Figure S1 The medianRank and Zsummary statistics of module preservation.

Notes: Module preservation was evaluated by medianRank and Zsummary statistics, which correlated to connectivity and density of networks. If Zsummary ⩾10, there is a strong evidence that the module is preserved. The module with a lower medianRank tends to exhibit stronger observed preservation than the module with a higher medianRank if both of them are preserved. (A) Module preservation analysis for the process of normal to dysplasia. (B) Module preservation analysis for the process of dysplasia to neoplasm. The pink module was found to be highly preserved in both processes of normal to dysplasia and dysplasia to OSCC.
Abbreviation: OSCC, oral squamous cell carcinoma.

Figure S2 (A) Volcano plot visualizing DEGs in GSE30784 (45 normal samples and 167 OSCC). The vertical lines demark the fold change values. The right vertical line corresponds to more than or equal to twofold up changes, and the left vertical line corresponds to twofold down changes, while the horizontal line marks a –log₁₀ adjusted P-value of 0.01. (B) Heat map hierarchical clustering reveals DEGs in OSCC groups compared with control groups.

Abbreviations: DEG, differently expressed genes; OSCC, oral squamous cell carcinoma.

Figure S3 Survival analysis of the 10 hub genes using GSE41613 for validation.

Table S1 GO enrichment analysis of DEGs

Category	Term	Count	%	P-value
Upregulated genes
BP	GO:0030198~ECM organization	60	8.39	6.57E–36
BP	GO:0007155~cell adhesion	78	10.91	9.15E–28
BP	GO:0030574~collagen catabolic process	32	4.48	8.50E–27
BP	GO:0006955~immune response	63	8.81	1.59E–19
BP	GO:0006954~inflammatory response	59	8.25	3.88E–19
BP	GO:0022617~ECM disassembly	25	3.50	8.43E–16
BP	GO:0035987~endodermal cell differentiation	16	2.24	8.80E–15
BP	GO:0009615~response to virus	28	3.92	1.36E–14
BP	GO:0060337~type I interferon signaling pathway	22	3.08	2.17E–14
BP	GO:0030199~collagen fibril organization	18	2.52	2.68E–14
CC	GO:0005576~extracellular region	175	24.48	4.88E–39
CC	GO:0031012~extracellular matrix	63	8.81	6.30E–29
CC	GO:0005615~extracellular space	140	19.58	1.26E–28
CC	GO:0005578~proteinaceous ECM	58	8.11	5.00E–27
CC	GO:0005581~collagen trimer	30	4.20	2.27E–19
MF	GO:0005201~ECM structural constituent	23	3.22	3.95E–15
MF	GO:0005518~collagen binding	21	2.94	5.06E–14
MF	GO:0008009~chemokine activity	18	2.52	1.92E–12
MF	GO:0005178~integrin binding	23	3.22	8.68E–11
MF	GO:0008201~heparin binding	28	3.92	1.25E–10
Downregulated genes
BP	GO:0030049~muscle filament sliding	12	1.76	3.40E–08
BP	GO:0006805~xenobiotic metabolic process	15	2.21	3.32E–07
BP	GO:0008544~epidermis development	15	2.21	9.90E–07
BP	GO:0006936~muscle contraction	16	2.35	3.36E–06
BP	GO:0055114~oxidation–reduction process	42	6.18	1.37E–05
BP	GO:0030855~epithelial cell differentiation	12	1.76	2.27E–05
BP	GO:0060048~cardiac muscle contraction	9	1.32	1.22E–04
BP	GO:0016338~calcium-independent cell–cell adhesion	6	0.88	5.83E–04
	via plasma membrane cell-adhesion molecules
BP	GO:0006941~striated muscle contraction	5	0.74	7.45E–04
BP	GO:0016266~O-glycan processing	9	1.32	9.28E–04
CC	GO:0070062~extracellular exosome	186	27.35	3.12E–20
CC	GO:0005615~extracellular space	89	13.09	2.98E–09
CC	GO:0016324~apical plasma membrane	29	4.26	9.39E–07
CC	GO:0005859~muscle myosin complex	7	1.03	9.48E–06
CC	GO:0030018~Z disc	16	2.35	1.25E–05
MF	GO:0003779~actin binding	32	4.71	8.33E–09
MF	GO:0008307~structural constituent of muscle	13	1.91	1.04E–08
MF	GO:0005198~structural molecule activity	24	3.53	1.51E–05
MF	GO:0016491~oxidoreductase activity	20	2.94	6.25E–05
MF	GO:0051371~muscle alpha-actinin binding	5	0.74	1.51E–04

Abbreviations: GO, Gene Ontology; DEG, differently expressed genes; BP, biological process; ECM, extracellular matrix; CC, cellular component; MF, molecular function.

Table S2 KEGG pathway enrichment analysis of turquoise module and brown module

Category	Term	Count	%	P-value
Upregulated genes
KEGG	hsa04512: ECM–receptor interaction	33	4.62	6.07E–21
KEGG	hsa04510: focal adhesion	46	6.43	9.77E–19
KEGG	hsa04151: PI3K–Akt signaling pathway	50	6.99	2.07E–12
KEGG	hsa05323: rheumatoid arthritis	24	3.36	9.34E–12
KEGG	hsa05146: amoebiasis	26	3.64	1.35E–11
KEGG	hsa04974: protein digestion and absorption	23	3.22	6.97E–11
KEGG	hsa04060: cytokine–cytokine receptor interaction	35	4.90	2.31E–09
KEGG	hsa05200: pathways in cancer	45	6.29	6.30E–08
KEGG	hsa05144: malaria	13	1.82	2.08E–06
KEGG	hsa04062: chemokine signaling pathway	26	3.64	2.09E–06
Downregulated genes
KEGG	hsa01100: metabolic pathways	68	10.00	1.80E–06
KEGG	hsa04530: tight junction	17	2.50	5.91E–06
KEGG	hsa04972: pancreatic secretion	13	1.91	3.28E–05
KEGG	hsa00982: drug metabolism – cytochrome P450	11	1.62	4.84E–05
KEGG	hsa00340: histidine metabolism	6	0.88	6.44E–04
KEGG	hsa00590: arachidonic acid metabolism	9	1.32	6.93E–04
KEGG	hsa00350: tyrosine metabolism	7	1.03	7.11E–04
KEGG	hsa05204: chemical carcinogenesis	10	1.47	8.98E–04
KEGG	hsa00980: metabolism of xenobiotics by cytochrome P450	9	1.32	0.002231
KEGG	hsa04670: leukocyte transendothelial migration	11	1.62	0.004088

Abbreviations: KEGG, Kyoto Encyclopedia of Genes and Genomes; ECM, extracellular matrix.