Identification and Verification of Five Potential Biomarkers Related to Skin and Thermal Injury Using Weighted Gene Co-Expression Network Analysis

Abstract

Background: Burn injury is a life-threatening disease that does not have ideal biomarkers. Therefore, this study first applied weighted gene co-expression network analysis (WGCNA) and differentially expressed gene (DEG) screening methods to identify pivotal genes and diagnostic biomarkers associated with the skin burn process. Methods: After obtaining transcriptomic datasets of burn patient skin and normal skin from Gene Expression Omnibus (GEO) and performing differential analysis and functional enrichment, WGCNA was used to identify hub gene modules associated with burn skin processes in the burn patient peripheral blood sample dataset and determine the correlation between modules and clinical features. Enrichment analysis was performed to identify the functions and pathways of key module genes. Differential analysis, WGCNA, protein-protein interaction analysis, and enrichment analysis were utilized to screen for hub genes. Hub genes were validated in two other GEO datasets, tested by immunohistochemistry for hub gene expression in burn patients, and receiver operating characteristic curve analysis was performed. Finally, we constructed the specific drug activity, transcription factors, and microRNA regulatory network of the five hub genes. Results: A total of 1,373 DEGs in GSE8056 were obtained, and the top 5 upregulated genes were S100A12, CXCL8, CXCL5, MMP3, and MMP1, whereas the top 5 downregulated genes were SCGB1D2, SCGB2A2, DCD, TSPAN8, and KRT25. DEGs were significantly enriched in the immunity, epidermal development, and skin development processes. In WGCNA, the yellow module was identified as the most closely associated module with tissue damage during the burn process, and the five hub genes (ANXA3, MCEMP1, MMP9, S100A12, and TCN1) were identified as the key genes for burn injury status, which consistently showed high expression in burn patient blood samples in the GSE37069 and GSE13902 datasets. Furthermore, we verified using immunohistochemistry that these five novel hub genes were also significantly elevated in burn patient skin. In addition, MCEMP1, MMP9, and S100A12 showed perfect diagnostic performance in the receiver operating characteristic analysis. Conclusion: In conclusion, we analyzed the changes in genetic processes in the skin during burns and used them to identify five potential novel diagnostic markers in blood samples from burn patients, which are important for burn patient diagnosis. In particular, MCEMP1, MMP9, and S100A12 are three key blood biomarkers that can be used to identify skin damage in burn patients.

Keywords: ROC; WGCNA; burn injury; peripheral blood; skin wound.

FIGURE 1

Study flowchart. (Study workflow. FC, fold change; GEO, Gene Expression Omnibus; GO, Gene Ontology; GS, gene significance; KEGG, Kyoto Encyclopedia of Genes and Genomes; MM, module membership; WGCNA, weighted gene co-expression network analysis).

FIGURE 2

Heatmap, volcano plot, and chromosome circos plot for differentially expressed genes identified in the GSE8056 dataset. The volcano map (A) indicates the difference of up-and downregulated genes, where red represents upregulation and blue represents downregulation. Heatmap (B) represents the differentially expressed genes expression patterns of differentially expressed genes in the upper burn patients versus the normal population group. Circos plot (C) shows the top 100 differentially expressed genes expression patterns and the distribution of the chromosomal location where they are located, with the outer circle representing the chromosome and the location of the gene in the chromosome, and the heatmap in the inner circle representing the expression of the top 100 differentially expressed genes (DEGs) in the burn dataset GSE8056. The top 5 upregulated differentially expressed genes (red) and top 5 downregulated differentially expressed genes (blue) according to |log2FC| values are connected in red and blue in the center of the circos plot, respectively.

FIGURE 3

GO and KEGG pathway analyses of DEGs identified in the GSE8056 dataset. (A) The histogram represents the GO pathway analysis of DEGs. The green bar represents BP (Biological process), blue represents CC (Cellular Component), and orange represents MF (Molecular Function). (B) Bubble diagram of the Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of differential genes.

FIGURE 4

Weighted gene co-expression network analysis (WGCNA) of burn-related key modules with PPI network analysis of key modules. (A) Sample clustering dendrogram of GSE19743 to detect outliers. (B) Analysis of the scale-free fit index (left) and the mean connectivity (right) for selecting various soft-thresholding powers (β). (C) Clustering dendrogram for genes in burn traits; each color below represents one co-expression gene module. (D) Heatmap depicting correlations between module and burn traits. (E) Scatter plot of the key module. Each point in the scatter plot represents one gene. (F) Hub genes in yellow module revealed by PPI using the cytoscape software.

FIGURE 5

The GO enrichment results of different modules in WGCNA are shown in the histogram, with different colors representing different categories of genes. Green represents BP, blue represents CC, and orange represents MF (BP: Biological process, CC: Cellular component, MF: Molecular function). (A) GO enrichment analysis results of the blue module. (B) GO enrichment analysis results of the brown module. (C) GO enrichment analysis results of the green module. (D) GO enrichment analysis results of the turquoise module. (E) GO enrichment analysis results of the yellow module.

FIGURE 6

Histogram of KEGG enrichment results of the modules with the count value and significant p-value. (A) Results of KEGG enrichment analysis for the blue module. (B) Results of KEGG enrichment analysis for the brown module. (C) Results of KEGG enrichment analysis for the green module. (D) Results of KEGG enrichment analysis for the turquoise module. (E) Results of KEGG enrichment analysis for the yellow module.

FIGURE 7

The ROC curve of five hub genes in GSE37069 and GSE139028. (A–E) The ROC curve of ANXA3, MCEMP1, MMP9, S100A12, and TCN1 in GSE37069. The x-axis shows specificity, and the y-axis shows sensitivity. ROC, receiver operating characteristic; AUC, area under the ROC curve.

FIGURE 8

Validation of expression levels of 5 hub genes in 4 burn-related GEO datasets. (A) The expression level of ANXA3, MCEMP1, MMP9, S100A12, and TCN1 between burned patients and normal patients in GSE8056. (B) The expression level of ANXA3, MCEMP1, MMP9, S100A12, and TCN1 between burned patients and normal patients in GSE19743. (C) The expression level of ANXA3, MCEMP1, MMP9, S100A12, and TCN1 between burned patients and normal patients in GSE37069. (D) The expression level of ANXA3, MCEMP1, MMP9, S100A12, and TCN1 between burned patients and normal patients from GSE139028. The red boxplot indicates the burn patient group, and the blue boxplot indicates the normal sample group. A t-test was performed to compare the means of the two groups (* represents p < 0.05, ** represents p < 0.01, *** represents p < 0.001, **** represents p < 0.0001, ns represents not significant).

FIGURE 9

Expression of the hub genes in burn injury skin tissues. (A) Representative images of immunohistochemical expression of ANXA3, MCEMP1, MMP9, S100A12, and TCN1 in deep second-degree burn skin tissues from second-degree burn patients (N = 13). (B) The histogram demonstrated that ANXA3, MCEMP1, S100A12, TCN1, and MMP9 showed positive expression in second-degree burn samples.

FIGURE 10

Construction of interaction network maps with transcription factors (A), miRNAs (B), and drug activity (C) for these five hub genes ANXA3, MCEMP1, MMP9, S100A12, and TCN1. (A) Interaction network diagram of the five hub genes and transcription factors, where red nodes represent key genes and blue nodes represent the transcription factors corresponding to the five hub genes. (B) Interaction network diagram of key genes and miRNAs, where red nodes represent key genes and yellow nodes represent miRNAs corresponding to the five hub genes. (C) Interaction network diagram of key genes and drug activities, where red nodes represent key genes and yellow nodes represent the drugs corresponding to the five hub genes.