Comprehensive Analysis of a Competing Endogenous RNA Co-Expression Network in Chronic Obstructive Pulmonary Disease

Introduction

Chronic obstructive pulmonary disease (COPD) is a complex respirator disorder featured by impaired expiratory airflow caused by airway and/or alveolar damage.^1,2 Statistics show that COPD is responsible for more than three million deaths worldwide each year.³ In the latest national COPD epidemiological survey, the total prevalence of COPD as measured by spirometry was 8.6%, which is equivalent to 99.9 million patients with COPD in China.⁴ Previous research has suggested that genetic susceptibility may influence COPD immune response and increase the risk of COPD under the influence of environmental exposures.^5–7 Further elucidating COPD pathogenesis and finding new biomarkers to improve treatment are important.

Non-coding RNAs (ncRNAs) are a class of transcripts that have been identified in the eukaryotic genome. Evidence suggests that ncRNAs are involved in numerous vital gene expression processes, such as chromatin remodelling/epigenetic memory, transcription, RNA splicing, editing, translation, and signal transduction.^8,9 lncRNAs and miRNAs have been implicated in COPD.¹⁰ The competitive endogenous RNA (ceRNA) hypothesis postulates that specific ncRNAs, for example, lncRNAs and circular RNAs (circRNAs), transcripts with shared microRNA (miRNA) binding sites competing for post-transcriptional control, may be involved in the downregulation of mRNA gene expression.¹¹ lncRNAs and their ceRNA-network have been demonstrated to be involved in the downregulation of mRNA gene expression in lung disorders.^12,13 SNHG5-miR-132-PTEN and CASC2-miR-18a-5p-IGF1 axis may play a potential role in COPD,^14,15 but ceRNA was still not explicated in COPD.

COPD induces lung tissue remodeling through sustained inflammatory responses of both congenital and adaptive immunity, in which neutrophils, macrophages, T lymphocytes and other inflammatory cells were involved.^16,17 Infection, inflammation, tobacco smoke, and immune response were able to provoke pulmonary autoimmunity and accelerate the pathogenesis of COPD.¹⁸ Yet, the precise immunological pathway involved remains indistinct.

This article aimed to explore an integrated transcriptome bioinformatics analysis based on mRNA and miRNA sequences to identify the disease-associated ceRNA-network and related genes. Considering the potential role of immune cells, we linked immune status with COPD and explored its value in the microenvironment of the disease.

Materials and Methods

The GSE38974 data set was used to obtain the RNA sequencing results generated from lung tissue samples. The mRNA expression profiles including COPD patients (N = 23) and healthy smokers (N = 9), based on the GPL4133 platform, were enrolled for analysis. The miRNA expression profiles including COPD patients (N = 19) and healthy smokers (N = 8), based on the GPL16384 platform, were enrolled. All patients and healthy smokers in this study were active smokers. Our research was performed according to the flowchart presented in Figure 1. The flowchart of analysis thoughts was shown in Figure S1.

Figure 1 Flow chart of the data analysis. Abbreviations: DEG, differentially expressed genes; DEmiRs, differentially expressed miRNAs.

Differential Expression Analysis

The GSE38974 data quality control was verified by principal component analysis (PCA). Differently expressed mRNAs (DEGs) and microRNAs (DEmiRs) were selected by difference analysis using the R package “limma” version 4.1.2. A significance threshold of p-value <0.05 and a fold change of | log2FC | >1 were adopted as the criteria for DEGs and DEmiRs. Furthermore, the “heatmap” and “ggplot2” packages were used to generate the heatmap, volcano plot, and box plot.

Functional Enrichment Analysis of DEGs

Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and Gene Set Enrichment Analysis (GSEA) pathway enrichment analysis were performed to determine the functional and molecular features of DEGs. The “ClusterProfiler” package of R software was used for enrichment analysis, in which preliminary analysis of the obtained DEGs was performed with systematic and comprehensive biological function notes attached.¹⁹ The GO term included biological processes (BP), cell components (CC), and molecular functions (MF). GO terms and KEGG pathways were deemed significant at a p-value of <0.05. For the visualization of the results, the R packages “enrichplot” and “GOplot” were utilized.

Construction of ceRNA-Network

Briefly, the construction of the ceRNA-network related to COPD can be divided into three stages. (a) Prediction of lncRNA–miRNA interactions: miRNA–lncRNA interactions can be predicted using the databases of starBase (v2.0),¹⁹ DIANA-LncBase (v3.0),²⁰ Encyclopedia of RNA Interactomes (ENCORI) database and lncACTdb;²¹ (b) Prediction of miRNA–mRNA interactions: miRNA–mRNA interactions were obtained from shared miRTarBase 2020²² and lncACTdb. miRTarBase predicted high-confidence functional miRNA–target interactions validated by reported assay and/or Western blot data. The miRNA–mRNA interactions were used to intersect with DEGs of GSE81292; (c) Construction of ceRNA-network: the predicted lncRNA–miRNA and miRNA–mRNA interactions were used to construct the lncRNA–miRNA–mRNA network. We used the R package “ggalluvial” to visualize the ceRNA-network by a Sankey plot.

Construction of a Weighted Gene Co-Expression Network

Weighted gene co-expression network analysis (WGCNA) is one of the most commonly used methods to construct modules and specific traits co-expression networks. The most variable upper 25% genes in GSE38974 mRNA and the whole gene of GSE38974 miRNA were selected for analysis.

First, we clustered the samples to evaluate whether there were any apparent outliers. Second, a weighted correlation network was calculated by constructing a pairwise correlation matrix between all gene pairs selected based on the variance. We chose the soft threshold power β values 5 for mRNAs and 4 for miRNA when the square value of the correlation coefficient reached 0.9 at the lowest power to build the scale-free topology model. Third, the Topological Overlap Measure was counted to reflect the overlap degree of shared neighbors, and we used hierarchical clustering and the dynamic tree cutting functions to perform the detection of module. Lastly, the module eigengene (ME) was defined as the first principal component of the gene module and can represent each module’s gene expression profiles.²³ Gene significance and module membership linked modules with medical features. The top 30 genes in the modules related to the traits were hub genes, and it was built by STRING with visualization by Cytoscape.

Assessing the Distribution of Immune Cell Subtypes

CIBERSORT is a powerful immune infiltration analysis tool based on the deconvolution integration method,²⁴ which is utilized to calculate the composition of varying immune cell species between COPD and control groups. This tool transforms a standardized gene expression matrix into the composition indices of 22 immune cell types. This algorithm was combined with the LM22 signature, and a sample of 1000 permutations with P < 0.05 was chosen for further analysis. The fractions of the 22 types of immune cells evaluated in each sample added up to 1. To display the ratio of different types of immune cells between COPD and control groups, a heat map was drawn via R software’s “pheatmap” package and a boxplot was generated using the R package “ggboxplot”. Additionally, the Wilcox rank-sum test was used to assess the association between the immune cells and the severity of Global Initiative for Chronic Obstructive Lung Disease (GOLD) stage.

Results

Differentially Expressed Genes

PCA was employed to assess the repeatability and quality of data within the groups. The repeatability and quality of data in GSE38974 was fine. In total, 622 DEGs (323 down-regulated and 299 up-regulated) were recognized in Figure 2A and B. 21 DEmiRNAs (14 down-regulated and 7 up-regulated) were recognized in Figure 2C and D. | log2FC| >1 and P-value <0.05 were used for the identification. DEGs functional enrichment analysis.

Figure 2 Differential expression analysis. (A) Volcano plots of DEGs; (B) Heatmap plots of DEGs; (C) Volcano plots of DEmiRs; (D) Heatmap plots of DEmiRs; (E) The top 15 KEGG pathways enrichment analysis of DEGs. The size of each circle represents the number of genes that are enriched in a particular pathway. The different color of each circle represents P-value. (F) The top 20 GO functional enrichment analyses of DEGs include biological process (BP), cellular component (CC), and molecular functions (MF) terms. The length of each bar represents the number of genes. The different color of each bar represents the analysis type. (G) GO terms Chord graph of DEGs. Abbreviations: DEG, differential expressed genes; DEmiRs, expressed miRNAs; KEGG, Kyoto Encyclopedia of Genes and Genomes; GO, gene ontology.

As shown in Figure 2E–G, the enrichment result of GO-BP terms included “positive regulation of RNA polymerase II promoter transcription”, “signal transduction”, “negative regulation of RNA polymerase II promoter transcription”, and “regulation of cell–cell adhesion”. The top enrichment result of GO-CC terms involved “extracellular region”, “extracellular space”, and “cytoplasm”. And “receptor-ligand activity”, “signaling receptor activator activity”, and “DNA-binding transcription activator activity, RNA polymerase II-specific” were the main enrichment results of GO-MF terms. KEGG pathway analysis indicated that all DEGs were primarily clustered in “cytokine-cytokine receptor interaction”, “lipid and atherosclerosis”, “tuberculosis”, and “transcriptional mis-regulation in cancer”. The most significant of these enriched GO terms and KEGG pathways indicated the functional interactions of DEGs. The results of up- and down-regulated DEGs GO and KEGG analyses were presented in Figure S2. The GSEA revealed that the two major pathways involved are the Cytokine-Cytokine Receptor Interaction and the JAK-STAT Signaling Pathway, which was shown in Figure S3. The top 20 pathways involved in GSEA KEGG enrichment analysis were listed in Table S1. The hub genes of the pathways were determined using the Cystoscope plug-in “cytohubba” and listed in Table S2.

Construction of a ceRNA-Network

We further short-listed 12 lncRNAs and 58 mRNAs by the intersections predicted by the abovesaid database. mRNAs or lncRNAs without the same prediction results were discarded (Figure 3A and B). Considering that there was a certain association between lncRNAs and mRNAs, those that are were not related will also be filtered by the lnACTdb dataset. In the issue, a co-expression network of COPD comprising 11 lncRNAs, 5 miRNAs, and 16 mRNAs was established and graphically depicted via a Sankey plot (Figure 3C).

Figure 3 Construction of the ceRNA Co-expression network. (A) Predicted mRNAs in Venn plot the intersection of DEGs, miRTarBase, and LncACTdb; (B) Predicted lncRNAs in Venn plot by the overlap of starBase, LncBase, ENCORI, and LncACTdb; (C) The Sankey plot of ceRNA network constructed by lncRNA–miRNA–mRNA.

WGCNA and Modules Construction

The top 25% of mRNAs with the most varying (N = 4937) were chosen for WGCNA network construction. If the genes showed similar expression patterns, they were divided into the same modules by hierarchical average linkage clustering. The selected soft threshold power was shown in Figure 4A for mRNAs. A dynamic tree-cutting algorithm and a merging threshold of 0.25 were utilized to identify eight modules (Figure 4B and C). Notably, the turquoise module was significantly correlated with COPD disease status (r = −0.88, P < 0.01). Figure 4D illustrated the significant association between the MM of the turquoise module and GS for COPD.

Figure 4 Construction of weighted gene co-expression network of mRNA. (A) The scale-free topology model fitting index for soft threshold powers (β) and the average connectivity for soft threshold powers. (B) Gene clustering tree dendrogram, with dissimilarity based on topological overlap in the mRNAs. Each color means one co-expressed module. (C) The heatmap of the 8 module which shows the relationship between module eigengenes and clinical traits. P-value was shown in each color cell. (D) Scatter plot of gene significance and module membership in the turquoise module.

Similarly, all miRNAs were used for WGCNA network construction. The soft-threshold power of 4 was selected, as depicted in Figure 5A. Subsequently, a cluster dendrogram plot of the gene network was created and displayed in Figure 5B. This enabled the identification of six co-expression modules. Module trait correlation analysis demonstrated that the green module (r = 0.97, P < 0.01) exhibited a statistically significant association with COPD (Figure 5C). Additionally, Figure 5D featured the significant correlation between MM and GS for COPD in the green module.

Figure 5 Construction of weighted gene co-expression network of miRNA. (A) Analysis of the scale-free topology model fitting index for soft threshold powers (β) and the average connectivity for soft threshold powers. (B) Gene clustering tree dendrogram, with the dissimilarity of topological overlap in the miRNAs. Each color means one co-expressed module. (C) The heatmap of 6 modules shows the relationship between model eigengenes and clinical traits. P-value is shown in each color cell. (D) Scatter plot of gene significance and module membership in green module.

Sub-Networks

The interaction network of the 10 mRNAs in the ceRNA-network was organized by STRING and visualization by Cytoscape (Figure 6A). The most significantly enriched functions of the 10 mRNAs were summarized in Figure 6B. The shared values of mRNAs in the ceRNA-network and hub-genes recognized by WGCNA including METTL7A, ZNF678, and PRRG4 were visualized in Figure 6C. Their binding miRNAs (hsa-miR-454-3p and has-miR-486-5p) were among the hub miRNAs identified by WGCNA and DEmiRNAs. Moreover, ZNF678, PRRG4, METTL7A, hsa-miR-454-3p and has-miR-486-5p were related to the GOLD stage and FVC group by linear Regression, respectively (Figures S4–7). The linear relationship between mRNA and miRNA was depicted in Figure S8. We filtered the ZNF678, PRRG4 and hsa-miR-454-3p and by Least absolute shrinkage and selection operator regression as well in Figure S9. Thus, one core sub-network: XIST, FGD5-AS1, KCNQ1OT1, HOXA11-AS, LINC00667, H19, PRKCQ-AS1, NUTM2A-AS1/hsa-miR-454-3p/ZNF678, PRRG4 was found (Figure 6D).

Figure 6 The core-ceRNA network of COPD. (A) PPI network of 10 mRNAs in the ceRNA network in Cytoscape software. (B) Heatmap of 10 mRNAs enrichment analysis results in ceRNA network. (C) Top 30 hub mRNAs in the weighted gene co-expression network turquoise module. The point where the mRNA in the ceRNA network overlap is marked in blue. (D) COPD’s core-ceRNA networks included XIST, FGD5-AS1, KCNQ1OT1, HOXA11-AS, LINC00667, H19, PRKCQ-AS1, NUTM2A-AS1/has-mir-454-3p/ZNF678, PRRG4. PPI, protein–protein interaction.

To further confirm the differential expression levels of hub genes between the control and COPD groups, GSE47460, containing 91 healthy samples and 145 COPD samples, and GSE151052, including 40 healthy samples and 77 COPD samples, were employed. Figure S10 shows that the expression levels of PRRG4 were significantly different between the two groups in the GSE47460 dataset. The expression level of ZNF678 was not measured in this dataset. In the GSE151052 dataset, PRRG4, and ZNF678 were differentially expressed between the control and COPD groups. These results verified that PRRG4 and ZNF678 were the hub genes of COPD.

Immunity Cells Analysis

The proportion of 22 immune cells in all profiles was estimated using the CIBERSORT algorithm, where the p-values obtained for all profiles were <0.05. The proportions of 22 immune cells for all profiles were presented in Figure 7A. Monocytes, T cells CD4 naïve, M2 macrophages, and M0 macrophages constituted the preponderance (Figure 7B and C). In addition, a comparison between COPD and controls revealed differences in terms of immune cell proportions (Figure 7D). The boxplot demonstrated that the COPD group had higher proportions of monocytes and M0 macrophages, as well as lower proportions of resting T cell CD8, NK cells activated, M2 macrophages, and mast cell resting in comparison to the control group. Additionally, the proportions of T cell CD8, NK cells activated, and macrophages M0 were correlated with the GOLD stage of COPD.

Figure 7 The profile of immune infiltration in COPD. (A) Boxplot of 22 immune cells composition in all groups. (B) Heatmap of the 22 immune cells between COPD and healthy smokers as controls. (C) The proportion of 22 immune cells in all groups. The column represents the samples, and the row represents the relative percentages of immune cells. (D) The boxplot shows the composition of different immune cell types between COPD and controls. Red represents COPD patients, while blue represents controls. *p<0.050, ** p<0.01, ***p<0.001.

Discussion

COPD is a complicated disease that damages both lung parenchyma and extrapulmonary systems.²⁵ Epigenetic mechanisms may play a role in the pathophysiology of COPD and may help to propose potential therapeutic targets.^26,27 Compared with proteome biomarkers such as high-degree complex and dynamic proteoforms, RNAs, including mRNAs, lncRNAs, and miRNAs in a transcriptome, are more straightforward and relatively stable, which are able to be the RNA biomarkers for effective COPD management.^28,29 In our research, expression profile data derived from the GEO dataset were compared between COPD patients and healthy smokers. To ensure reliable results, hub genes in COPD were identified using the WGCNA to identify connected nodes. Subsequently, a ceRNA-network was constructed, comprising 11 lncRNAs, 5 miRNAs, and 16 mRNAs. One core sub-network: XIST, FGD5-AS1, KCNQ1OT1, HOXA11-AS, LINC00667, H19, PRKCQ-AS1, NUTM2A-AS1/has-mir-454-3p/ZNF678, PRRG4 was obtained.

The COPD patients’ gene expression profiles were previously studied to explore the gene markers from public databases.^30–32 ceRNA-network via DEGs analysis in primally of human small airway was constructed.³³ In addition, DEGs analysis based on multiple gene expression dataset was also used to construct the ceRNA-network by mRNA expression data.³² However, the confounding factors may decrease the sensitivity and the specificity of the DEGs as the biomarkers for COPD. In comparison to other bioinformatics approaches, WGCNA results are more reliable and have greater clinical value due to the focus on correlations between clinical traits and co-expression modules.²³ Genes clustered in the same module are presumed to be functionally related, thus making the outcome more dependable when the common value of DEGs and hub genes obtained through WGCNA is applied. By leveraging biologically related modules and hub genes, there is considerable potential for them to serve as predictive biomarkers for diagnosis or treatment.

The ceRNA-network hypothesis proposes a regulatory mechanism mediated by lncRNAs. In this theory, lncRNAs are used as an endogenous miRNA sponge to regulate mRNAs expression, which has enhanced our comprehension of the transcriptional-level mechanisms of some illnesses but still has disputes in COPD.³⁴ We selected the interaction of miRNA–mRNA by using LncACTdb and miRTarBase database. Owing to the dearth of lncRNA sequencing data of COPD, multiple databases, such as ENCORI, starBase, LncACTdb, and LncBase, were employed to predict miRNA–lncRNA interactions. These databases are built on experimentally supported evidence or statistically predictions, and thus can be deemed reliable. Eventually, we obtained one core ceRNA-network: XIST, FGD5-AS1, KCNQ1OT1, HOXA11-AS, LINC00667, H19, PRKCQ-AS1, NUTM2A-AS1/has-mir-454-3p/ZNF678, PRRG4. The memberships of the ceRNA-network were potentially involved in various inflammation diseases.^35,36 XIST was also proved to exacerbate rheumatoid arthritis and systemic lupus erythematosus.^37,38 The function of lncRNAs as ceRNAs network was also verified to regulate miRNA–gene interactions in COPD.^14,15

The impact of miR-454-3p on inflammation has been widely studied. Notably, lower expression of miR-454-3p was associated with decreased inflammation and regression of cirrhosis during treatment of hepatitis B virus infection.³⁹ Inhibition of miR-454-3p is thought to induce M2 macrophage polarization in glioma.⁴⁰ Meanwhile, the miR-454-3p suppression attenuated the viability of human aortic endothelial cells, while overexpression of miR-454-3p inhibited the apoptosis of HAEC-induced ox-LDL, thereby affecting the apoptosis of endothelial cells.⁴¹ Transmembrane gamma-carboxyglutamic acid protein 4 (PRRG 4), a protein-coding gene has been proved to promote the progression of breast cancer, and enhance metastasis.⁴² It can also promote autism symptoms of WAGR by perturbing these mechanisms in the developing human brain.⁴³ A recent study also verified that PRRG4 was the main methylation marker of COPD,⁴⁴ which is verified by our study as well. The potential mechanism of ZNF678 in COPD remains unclear. The susceptibility and severity of COPD vary between different genders.⁴⁵ In this study, we conducted an analysis of the expression levels of hub mRNA and miRNA between male and female COPD patients. Our analysis revealed differential expression of has-miR-10a-5p and CCDC137 in male and female COPD patients (Figures S11 and 12).

This study revealed that, according to the results of the CIBERSORT algorithm, the lung tissues of COPD patients had higher proportions of monocytes, M0 macrophages, and resting T cell CD8 than healthy smokers, whereas the proportions of NK cells activated, M2 macrophages, and resting mast cell tended to be lower. And T cell CD8, NK cells activated, and macrophages M0 were associated with the severity of COPD. Recently, therapeutic medication that inhibited the inflammatory cascade, chemotaxis, or inflammatory cell activation plays a crucial role to delay the progression of airway remodeling in COPD and to promote the repair of the damaged airway.⁴⁶ Our data suggested that the differential immune cells in COPD in comparison to controls may provide new targets for the treatment of COPD.

This study was limited by several factors. Firstly, the data for analysis was sourced from GEO online data, thus requiring further experiments to verify the results and assess their underlying medical value. Secondly, only miRNA and mRNA expression profiles of COPD patients were identified, suggesting the need for further research to explore the underlying mechanisms of the disease. Further lncRNA is warranted to verify our results. Thirdly, we did not confirm the interactions between miRNA and mRNA in the ceRNA network. Finally, the potential diagnostic and therapeutic roles of the ceRNA-network in COPD were not explored.

In future research endeavors, we intend to investigate the expression levels of these markers in both in vivo and in vitro experiments, as well as their potential involvement in the pathogenesis of COPD. Furthermore, we may focus on the ceRNA-network member as the diagnostic markers and as therapeutic targets as well. The data have provided clues for further research on the regulation of the potential miRNAs and mRNAs treatment targets in COPD.

Conclusion

Overall, this study centered on understanding the regulatory role of the ceRNA-network in COPD and identified reliable core-ceRNA interaction pairs. The results of the immune infiltration analysis further suggest that immune cells in the lung may indicate the progression of the disease and be involved in the progression of COPD. The study may provide new evidence for the pathogenesis of COPD, although the exact molecular mechanism and regulatory pathways of COPD still need to be further explored.