Back to Journals » International Journal of Chronic Obstructive Pulmonary Disease » Volume 16

Identifying miRNA Modules and Related Pathways of Chronic Obstructive Pulmonary Disease Associated Emphysema by Weighted Gene Co-Expression Network Analysis

Authors An J, Yang T, Dong J, Liao Z, Wan C, Shen Y, Chen L 

Received 24 June 2021

Accepted for publication 25 October 2021

Published 15 November 2021 Volume 2021:16 Pages 3119—3130


Checked for plagiarism Yes

Review by Single anonymous peer review

Peer reviewer comments 2

Editor who approved publication: Professor Zhang

Jing An, Ting Yang, Jiajia Dong, Zenglin Liao, Chun Wan, Yongchun Shen, Lei Chen

Department of Respiratory and Critical Care Medicine, Division of Pulmonary Diseases, State Key Laboratory of Biotherapy, West China Hospital, Sichuan University, Chengdu, Sichuan, 610041, People’s Republic of China

Correspondence: Yongchun Shen; Lei Chen Email [email protected]; [email protected]

Background: Chronic obstructive pulmonary disease (COPD) is a heterogeneous chronic inflammatory disease characterized by progressive airflow limitation that causes high morbidity and mortality. MicroRNA, a short-chain noncoding RNA, regulates gene expression at the transcriptional level. microRNA modules with a role in the pathogenesis of COPD may serve as COPD biomarkers.
Methods: We downloaded the GSE33336 microarray data set from the Gene Expression Omnibus (GEO) database, the data are derived from 29 lung samples of patients with emphysema undergoing curative resection for lung cancer. We used weighted gene co-expression network analysis (WGCNA) to construct co-expression modules and detect trait-related microRNA modules. We used the Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis to predict the biological function of the interest modules, and we screened out candidate hub microRNAs based on their module membership (MM) value and top proteins on the results of the protein–protein interaction (PPI) network.
Results: Three microRNA modules (royal blue, light yellow and grey60) were highly associated with COPD. Axon guidance, proteoglycans in cancer and mitogen-activated protein kinases (MAPK) signaling pathway were common pathways in these three modules. Keratin18 (KRT18) was the top protein in our study. miR-452, miR-149, miR-133a, miR-181a and miR-421 in hub microRNAs may play a role in COPD.
Conclusion: These findings provide evidence for the role of miRNAs in COPD and identify biomarker candidates.

Keywords: chronic obstructive pulmonary disease, weighted gene co-expression network analysis, microRNA


Chronic obstructive pulmonary disease (COPD) is a major health burdens worldwide causing high morbidity and mortality.1 The disease is predicted to become one of the top four causes of death globally by 2030.2 COPD is a chronic inflammatory disease characterized by progressive airflow limitation involving peripheral loss of lung parenchymal tissue (emphysema) and mucus hypersecretion with chronic bronchitis.1 These pathophysiological changes result in a progressive decline of forced expiratory volume in the first second, after full inspiration (FEV1), inadequate lung emptying on expiration, and subsequent static and dynamic hyperinflation.3 COPD is a heterogeneous and complex disease caused by diverse cellular and pathophysiological changes under distinct genetic backgrounds.4 Studies have associated single genes, such as the gene encoding matrix metalloproteinase 12 (MMP-12) and glutathione S-transferase, to a decline in lung function5 and an increased COPD risk.6 Genome-wide association studies have linked genetic loci (including FAM13A, markers near the alpha-nicotinic acetylcholine receptor, hedgehog interacting protein (HHIP), and others) with COPD7–11 These gene expression and polymorphism studies have shown that COPD is a polygenic disease involving multiple pathogenic signaling pathways.

MicroRNAs (miRNAs) are a class of small, single‑stranded noncoding RNAs. By binding to the 3′terminal noncoding region of target mRNAs, miRNAs regulate gene expression at the transcriptional level in the form of complete complementation or incomplete complementation.12,13 miRNAs affect biological processes, including differentiation, proliferation, and apoptosis, through their regulation of gene expression under normal and pathological conditions.14,15 Studies have suggested pathophysiological roles for miRNAs in COPD: For example, Cao et al detected increased expression levels of miR-183, miR-200b, and miR-200c in the blood of patients with COPD16 and, Donaldson et al found that plasma levels of muscle-specific miR-499, miR-133, and miR-206 were elevated in patients with COPD.17 By using differential gene expression analysis, these studies focused on one or several miRNAs and only revealed simple up- and down-regulation, which cannot build the relationships between miRNAs and disease.

The weighted gene co-expression network analysis (WGCNA) is a systems biology method for describing the correlation patterns among genes across microarray samples.18 The method is based on the theory that genes with closely functional linkages or involved in similar pathways may have similar expression profiles.19,20 By using information from thousands of genes with greatest changes, WGCNA identifies interest gene clusters (which we call gene modules). Next, WGCNA identifies associations between modules and external clinical traits, providing an effective way to explore the mechanisms behind certain traits.18 WGCNA has been widely used in cancer research.21–23 Qin et al used WGCNA to investigate and contrast the molecular processes differing between the bronchiolitis and emphysema phenotypes of COPD.24

The aim of this study was to use the WGCNA method to construct co-expression modules for the miRNA expression data from human biological specimens and to identify correlations between different modules and clinical COPD traits. This methodology provides an effective way to explore the role of miRNAs from different perspectives and identifies promising diagnostic biomarkers and molecular COPD pathophysiological pathways.


Microarray Data Used

We downloaded the miRNA expression profile in the GSE33336 dataset from the Gene Expression Omnibus (GEO) public functional genomics data repository of the National Center of Biotechnology Information (, platform GPL6955. To obtain this dataset, researchers set up Agilent human miRNA oligo arrays (8x15K array, Part No. G4470A, Agilent Technologies) based on miRBase V9.1 on lung tissues from 29 emphysematous lungs obtained from patients undergoing curative resection for lung cancer.25 Briefly, smokers with at least a 20-pack-year history who quit smoking at least 10 months prior to surgery were considered. All patients had FEV1/ (forced vital capacity) FVC ratio <0.70, indicating the presence of COPD. Thus, 29 former smokers with COPD were selected and classified according to gas transfer measurement results (single breath carbon monoxide diffusion coefficient, KCO) as having “mild” emphysema (>75% predicted KCO) or “moderate” emphysema (40–75% predicted KCO). Specimens from patients on inhaled or oral steroids, with lung pathologies that might confound spirometry measurements, or with alpha 1 antitrypsin (α1AT) enzyme deficiency were excluded. Detailed information can be found in other publications.25,26

Identification of Differentially Expressed miRNAs (DEMs)

For the original array data, the context correction and normalization were performed using a robust multiarray average process. The data were subsequently converted into expressive measures using the R affy package. DEMs were identified using the Limma package. The cutoff values for DEM screening using the Benjamini & Hochberg procedure were a p ≤0.05 and an absolute Log2Ratio ≥1.


The WGCNA package is a comprehensive collection of R functions for weighted correlation network analysis.18,27 Datasets were constructed using a previously described method.18 First, we constructed a signed weighted correlation network by creating a matrix of pairwise correlations between all pairs of miRNAs chosen according to their variance. The resulting Pearson correlation matrix was transformed into a matrix of connection strengths (eg, an adjacency matrix). Next, we applied the R pickSoftThreshold function that performs an analysis of network topology and chose a proper soft-thresholding power of 4.28 The topological overlap was calculated to measure the network interconnectedness. Average linkage hierarchical clustering was used to group genes based on the topological overlap dissimilarity measure (1-topological overlap) of their network connection strengths. Using a dynamic tree-cutting algorithm and merging the threshold function at 0.25, we identified 25 modules in the datasets. We generated heatmap plots to visualize the co-expression module structure of topological overlap in the gene network. The associations among modules were summarized by a hierarchical clustering dendrogram of the eigengenes and by a heatmap plot of the corresponding eigengene network.

Association of Modules to External Traits

To identify modules that were significantly associated with the traits of gender, age, COPD status (mild or moderate), predicted FEV1, predicted KCO, and smoking pack years we correlated the module eigengenes (MEs; ie, the first principal component of a module)29 with external traits using Pearson’s correlation coefficients, modules with p-values <0.05 were identified as trait-related modules.

Identification of Hub miRNAs and Construction of Protein–Protein Interaction

Module membership (MM) represents the correlation between the gene expression profile and the module eigengene. We considered the top 40 miRNAs with the highest module membership values in the modules as hub miRNAs and visualized them with Cytoscape. We obtained the target genes from the database “microWalk” ( These target genes have been validated by existing experimental procedures. The target genes of hub microRNAs were used to construct a protein–protein interaction (PPI) network in the STRING database ( to identify associated functional proteins and important microRNAs.

Functional Annotation and Enrichment Analysis

We performed functional annotations with the ClusterProfiler based on the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases to assign biological functions to the miRNAs in these modules. Target genes for the miRNAs in each module were subjected to this analysis. We visualized the top 20 most enriched function terms (ordered by their p-values).


Identification of DEMs in COPD

The miRNA-sequencing data from the GEO dataset were subjected to differential expression analysis based on a threshold of p ≤0.05 and an absolute Log2Ratio ≥1. We found 422 DEMs that were significantly downregulated and 82 DEMs that were significantly upregulated in the samples from patients with moderate COPD samples when compared with the samples from patients with mild COPD (Table S2). The Additional file 1 (Figure S1) displays volcano plots of the miRNAs.

Identification of miRNA Coexpression Networks and Modules

We selected 9800 probes for network construction. miRNAs exhibiting similar patterns of expression were grouped into modules via hierarchical average linkage clustering. The sample dendrogram and trait heatmap grouped the selected samples into different clusters and provided the clinical trait data distribution map (Figure 1). Network topology was analyzed using various soft threshold powers at 4 (Additional file 2: Figure S2), according to the method by Zhang and Horvath.30 After using a dynamic tree cutting algorithm (Figure 1), we obtained 25 distinct co-expression modules; 43 uncorrelated miRNAs were assigned into a grey module ignored in the following study. We used the Pearson correlation coefficient to analyze the interactions of these co-expression modules (Figure 1). The darker background indicates a higher module correlation.

Figure 1 Sample clustering and module detection.

Notes: (A) Sample cluster tree and trait indicator. The leaves of the tree correspond to the sample. The red color represents female or male, and COPD status. The color intensity was proportional to older age, higher FEV1, higher KCO and higher smoke-pack-year. (B) Clustering dendrogram of miRNAs, with dissimilarity based on topological overlap, together with assigned module colors. (C) Visualizing the miRNA network using a heatmap plot. The heatmap depicts the Topological Overlap Matrix (TOM) among all miRNAs in the analysis, light color represents low overlap and progressively darker red color represents higher overlap. Blocks of darker colors along the diagonal are the modules. The miRNA dendrogram and module assignment are also shown along the left side and the top.

Association of miRNA Co-Expression Modules with Clinical Traits

We correlated the 25 modules with traits of interest and searched for the most significant associations to assess the physiological significance of the modules. In the heatmap of module-trait correlation (Figure 2), rows indicate modules and columns indicate external traits. miRNAs clustered in the royal blue module had the strongest negative correlation with the predicted FEV1 (correlation, −0.41; p-value, 0.02), miRNAs clustered in the light yellow module had the strongest positive correlation with the moderate COPD status (correlation, 0.46; p-value, 0.008), miRNAs clustered in the grey60 module had the strongest negative correlation with the smoking pack years (correlation, −0.61; p-value, 3e-04). Thus, we considered mainly modules with strong correlation to interesting traits to ensure biological significance. The royal blue, light yellow, and grey60 modules were analyzed further: The additional file 6 (Table S1) shows the miRNAs in these three modules associated to clinical traits of interest.

Figure 2 Module-trait associations. Each row corresponds to a module eigengene, and each column to a trait. Each cell contains the corresponding correlation and p-value. The table is color-coded by correlation according to the color bar.

Functional Annotation and Enrichment Analysis

We performed functional annotation and enrichment analyses to explore the biological functions associated with the genes clustered in each module of interest (royal blue, light yellow, and grey60). The GO database analysis resulted in three main annotated categories (biological processes, cellular components, and molecular functions). For each module and each function database, we generated bar plots to visualize the top 20 most significantly enriched terms. The additional files 3–5 (Figures S3S5) show the summarized results of the GO enrichment analysis with important findings (p-value < 0.05).

For the royal blue module, the top five significantly enriched biological processes pathways included intracellular signal transduction, intracellular protein transport, neutrophil degranulation, positive regulation of GTPase activity, and endocytosis. The top five significantly enriched cellular component pathways were glutamatergic synapse, neuron projection, membrane raft, neuronal cell body, and axon. Moreover, the top five significantly enriched molecular function component pathways were DNA-binding transcription activator activity (RNA polymerase II-specific), protein kinase binding, protein serine/threonine kinase activity, transcription regulatory region DNA binding, and protein domain specific binding. As shown in Figure 3, the KEGG annotation results showed axon guidance, proteoglycans in cancer, MAPK signaling pathway, hepatocellular carcinoma, and sphingolipid signaling pathway as the top five significantly enriched pathways.

Figure 3 Kyoto Encyclopedia of Genes and Genomes pathways enrichment. (A) royal blue module (B) light yellow module (C) grey60 module.

For the light yellow module, the top five significantly enriched biological processes pathways were intracellular signal transduction, intracellular protein transport, positive regulation of GTPase activity, brain development, and the Wnt signaling pathway. The top five significantly enriched cellular component pathways were glutamatergic synapse, neuronal cell body, neuron projection, membrane raft, and dendrites. Further, the top five significantly enriched molecular function component pathways were DNA-binding transcription activator activity (RNA polymerase II-specific), GTPase activator activity, Rab GTPase binding, protein serine/threonine kinase activity, and protein kinase binding. As shown in Figure 3, the KEGG annotation results showed MAPK signaling pathway, proteoglycans in cancer, axon guidance, hippo signaling pathway, and inositol phosphate metabolism as the top five significantly enriched pathways.

For the grey60 module, the top five significantly enriched biological processes pathways were intracellular signal transduction, positive regulation of GTPase activity, intracellular protein transport, regulation of small GTPase mediated signal transduction, and the Wnt signaling pathway. The top five significantly enriched cellular component pathways were glutamatergic synapse, neuron projection, neuronal cell body, early endosome membrane, and axon. In addition, the top five significantly enriched molecular function component pathways were GTPase activator activity, DNA-binding transcription activator activity, RNA polymerase II-specific, protein serine/threonine kinase activity, DNA-binding transcription repressor activity, RNA polymerase II-specific, and transcription regulatory region DNA binding. As shown in Figure 3, the KEGG annotation results showed proteoglycans in cancer, axon guidance, endocytosis, signaling pathways regulating pluripotency of stem cells, and MAPK signaling pathway as the top five significantly enriched pathways.

Identification of Hub miRNAs and PPI

Module memberships (MMs) represent each microRNA’s Pearson correlation coefficient with its corresponding eigengene module. We considered the top 40 miRNAs with the highest module membership values as hub miRNAs and used Cytoscape to visualize them (Figure 4 and Table S3). The hub microRNAs miR-452, miR-149, miR-133a, miR-181a, and miR-421 may play a role in emphysema-associated COPD. We used the target genes of these hub microRNAs for protein–protein interaction analysis to get a systematic perspective of the molecular mechanisms of this disease. In PPI network, the degree means the connection between proteins. The higher the degree, the stronger the connection is. Protein with highest degree is the top protein. As shown in Figure 5, the top three proteins in the royal blue module were KRT18, insulin-like growth factors-1 (IGF1), and KRT84 (their degrees were 26, 21, and 19, respectively). The top three proteins in the light yellow module were KRT71, KRT84, and KRT2 (their degrees were 11, 10, and 10, respectively). The top three proteins in the grey60 module were KRT18, KRT2, and KRT76 (their degrees were 21, 17, and 17, respectively). miR-520b targets the KRT18 protein in both the royal blue and the grey60 modules.

Figure 4 Visualization of the network connections among the most connected miRNAs in the royal blue module (A) light yellow module (B) and grey60 module (C) generated by the Cytoscape software. Edge weight indicates the TOM similarity between two nodes.

Figure 5 The protein–protein interaction (PPI) analysis of target gene of hub microRNAs in royal blue module (A) light yellow module (B) and grey60 module (C). Degrees describe the importance of the protein nodes.


For this study, we applied a systems biology approach (WGCNA) to identify miRNA modules and associated pathways of emphysema-associated COPD for the first time. We identified 25 modules by WGCNA analysis and found that the royal blue, light yellow, and grey60 modules were significantly correlated to our traits of interest (predicted FEV1, moderate COPD status, and smoking pack years). These three modules were analyzed further.

Axon guidance, proteoglycans in cancer, and the MAPK signaling pathway were common processes in the top 20 KEGG pathways of the three modules. Axon guidance plays an important role during the development of neural circuits. Several receptors and pathways have also been shown to play roles in axon guidance.31 Slit guidance ligand (SLIT) and roundabout (ROBO) were shown to be involved in neural development by inhibiting the migration of axons.32 Subsequently, SLIT was shown to be a potential anti-inflammatory molecule.33 Lin et al used an integrated bioinformatics analysis and found that the expressions of SLIT2 and ROBO2 genes are downregulated in patients with COPD and that the expressions are significantly negatively correlated with the COPD stages.34 The Wnt pathway is an evolutionarily conserved family that plays a crucial role in axon and dendrite development.35 Wnt signaling modulates several cellular processes involved in COPD, including oxidative stress, apoptosis/proliferation, inflammation, mucus hypersecretion, protease/antiprotease imbalance, autophagy, senescence, metabolism reprogramming, mitochondrial dysfunction and stem-cell/progenitor-cell renewal.36

Proteoglycans are a family of charged molecules that contain a core protein and at least one glycosaminoglycan. They are primarily found on the cell surface and in the extracellular matrix, where they serve as structural components. They are thought to influence a myriad of cellular processes, including those linked with cancer development.37 Proteoglycans are found in the extracellular matrix, in plasma membrane of cells, and as intracellular structures in the lungs.38 Changes in the proteoglycan composition in the lungs have been reported in human lung diseases, such as idiopathic pulmonary fibrosis, acute respiratory distress syndrome (ARDS), COPD, lymphangioleiomyomatosis, and asthma.39 Studies have shown that proteoglycans and their associated glycosaminoglycans play important roles modulating pulmonary inflammation.40–42 Takahashi et al found that proteoglycans maintain lung stability in an elastase-treated mouse model of emphysema, suggesting that the loss of proteoglycans in human emphysema contributes to disease progression.43

Mitogen-activated protein kinases (MAPKs) are serine-threonine protein kinases including c-Jun NH2-terminal kinase (JNK), p38 MAPK, and extracellular signal- regulated kinase (ERK). MAPKs regulate several cellular activities such as proliferation, differentiation, apoptosis, survival, inflammation, and innate immunity.44 MAPK signaling pathways have also been found to play key roles in the pathogenesis of many diseases.45 The JNK and p38 MAPK signaling pathways are activated by various types of cellular stress such as oxidative stress, bacterial lipopolysaccharide, and proinflammatory cytokines (tumor necrosis factor-α and interleukin 1β);44 these factors also play a role in the pathogenesis of COPD.4 Patients with COPD have increased p38 MAPK activation in macrophages and in other cells within the alveolar wall. Studies have shown that p38 inhibitors exhibit anti- inflammation effects both in vitro and in vivo.46 However, many p38 inhibitors have failed in clinical trials due to their unacceptable safety profiles.47

In addition to KEGG enrichment analysis of the key modules in our study, we used the STRING database to generate a protein–protein interaction (PPI) network of some hub genes in key modules. The PPI network analysis results showed that KRT18 had the highest degree value in all the module networks. Cytokeratins are major structural proteins of epithelia that are essential for tissue integrity.48 Cytokeratin 18 (CK18, also known as KRT18) is an acidic cytokeratin member of the protein family of cytoskeletal intermediate filaments.49 It is widely expressed in normal epithelia, including in pulmonary bronchial and alveolar epithelial cells.50,51 CK18 is involved in important signaling pathways including apoptosis, cell cycle, and cancer progression.49 During apoptosis, specific caspase-cleaved CK18 (ccCK18) fragments are generated. Hacker and Lichtenauer et al demonstrated that elevated levels of ccCK-18 in the serum of COPD patients correlate with disease severity, and that hypoxic conditions may trigger the release of the ccCK-18 fragments.52,53

Studies have shown that IGF-1 signaling plays an important role during lung development and diseases (such as congenital disorders, cancers, inflammation, and fibrosis).54 IGF-1 signaling is involved in the pathogenesis of COPD by promoting inflammation, muscle dysfunction metabolism, and the aging process.54

In this study, we identified several hub miRNAs that may play roles in the development of emphysema-associated COPD. miR-452, miR-149, miR-133a, miR-181a, and miR-421 had been found to be abnormally expressed in patients with emphysema,55 pulmonary fibrosis,56 pulmonary arterial hypertension,57 and lung cancer.58 The expression of miR-452 in alveolar macrophages is reduced in smokers; miR-452 targets matrix metalloproteinase-12 (MMP-12), which is associated with the development of emphysema.55 Shen et al reported that treatment with cigarette smoke extract (CSE) decreases the level of miR-149-3p in a murine monocytic cell line (THP-1): Transfection of miR-149-3p mimic into CSE-induced THP-1 cells decreased the levels of the TLR-4 and NF-κB p65 proteins, which play an important role in COPD.59 Transforming growth factor (TGF)-beta1-induced miR-133a could inhibit myofibroblast differentiation and pulmonary fibrosis.56 Moreover, myofibroblast differentiation has been observed in patients with COPD.60 Liu et al found that the expression of miR-133 was lower in patients with acute COPD exacerbations than in stable patients and healthy controls.61 Hypercapnia, a common occurrence in patients with COPD, was found to increase the airway smooth muscle contractility via caspase-7-mediated miR-133a-RhoA signals.62 miR-181a was found to target endocan and thereby ameliorate the inflammatory response in monocrotaline-induced pulmonary arterial hypertension.63 Studies have found that endocan is associated with disease severity, exacerbation, and lung function decline in patients with COPD.63–65 Yuan et al demonstrated that miR-421 induces apoptosis in the lung and increases the inflammatory response in a bronchopulmonary dysplasia mouse model.66 Knock-down of miR-421 causes levels of KEAP to rise, inhibiting Nrf2-dependent antioxidant expression, and ultimately increasing the intracellular ROS levels in non-small cell lung cancer.58 Nrf2 and oxidative stress play important roles in COPD pathogenesis.67 The role of miR-421 in COPD needs further exploration; but we hypothesize that the profiles of miR-452, miR-149, miR-133a, miR-181a, and miR-421 may be associated with the inflammation, myofibroblast differentiation, oxidative stress, and apoptosis evident in COPD. This hypothesis needs be confirmed by others.

We identified a few miRNAs related to in interesting clinical traits. Jing et al found that miR-558 could reduce the damage of HBE cells exposed to cigarette smoke extract (CSE) by targeting TNF receptor superfamily member 1A (TNFRSF1A) and inactivating TAK1/MAPK/NF-κB pathway.68 miR-558 is demonstrated to participate in the lung cancer cell growth, tumorigenesis, and invasion.69,70 Li et al demonstrated that miR −587 acts as an oncogene in non-small-cell lung carcinoma.71 Several studies showed that miR-621, miR −639, miR −566, miR −769, miR −597 play a role in proliferation, metastasis, migration, and invasion in a variety of cancer cells.72–76

We are aware of the limitations in this study. First, normal lung tissues were not included in the study and this omission may have caused biases. Second, we analyzed a single platform dataset; therefore, the relatively small sample size may not fully represent all patients with COPD. Third, all the COPD lung samples were from patients who also had lung cancer, and the presence of the cancer may have affected our results. Finally, we did not provide biological experimental verification of the hub miRNAs and their associations with genes. The results need to be verified by functional experiments as a next step of the research.


In this study, we identified three microRNA modules related to emphysema-associated COPD through WGCNA. Axon guidance, proteoglycans in cancer, and the MAPK signaling pathway may play a role in the pathological process of emphysema-associated COPD. KRT18, miR-133a, miR-181a, and miR-421 may be associated with emphysema-associated COPD according to our combined KEGG enrichment and PPI network analyses. These findings should be verified and extended in future experimental work.


COPD, chronic obstructive pulmonary disease; GEO, Gene Expression Omnibus; WGCNA, weighted gene co-expression network analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes enrichment analysis; MM, module membership; PPI, protein–protein interaction; FEV1, forced expiratory volume in 1 s; MMP-12, matrix metalloproteinase 12; HHIP, hedgehog interacting protein; FVC, forced vital capacity; KCO, single breath carbon monoxide diffusion coefficient; α1AT, alpha 1 antitrypsin; DEMs, differentially expressed miRNAs; MEs, module eigengenes; MAPK, mitogen-activated protein kinase; SLIT, Slit guidance ligand; ROBO, roundabout; ARDS, acute respiratory distress syndrome; JNK, c-Jun NH2-terminal kinase; ERK, extracellular signal- regulated kinase; KRT, keratin; IGF-1, insulin-like growth factors-1; TGF, Transforming growth factor; TLR-4, Toll-like receptor 4; KEAP, Kelch-like ECH-associated protein; Nrf2, nuclear factor E2-related factor 2; ROS, reactive oxygen species.

Data Sharing Statement

The datasets generated during and/or analyzed during the current study are available in the Gene Expression Omnibus (GEO) datasets (


This study was supported by the National Natural Science Foundation of China (31871157, 81830001, 81500010) and 1·3·5 project for disciplines of excellence, West China Hospital, Sichuan University (ZYGD18006, ZYJC18012).


The authors declare that they have no competing interests.


1. Rabe KF, Watz H. Chronic obstructive pulmonary disease. Lancet. 2017;389:1931‐1940. doi:10.1016/S0140-6736(17)31222-9

2. McLean S, Hoogendoorn M, Hoogenveen RT, et al. Projecting the COPD population and costs in England and Scotland: 2011 to 2030. Sci Rep. 2016;6:31893. doi:10.1038/srep31893

3. Regan EA, Lynch DA, Curran-Everett D, et al. Clinical and radiologic disease in smokers with normal spirometry. JAMA Intern Med. 2015;175:1539–1549. doi:10.1001/jamainternmed.2015.2735

4. Celli BR, Decramer M, Wedzicha JA, et al. An official American Thoracic Society/European Respiratory Society statement: research questions in COPD. Eur Respir Rev. 2015;24:159–172. doi:10.1183/16000617.00000315

5. McCloskey SC, Patel BD, Hinchliffe SJ, et al. Siblings of patients with severe chronic obstructive pulmonary disease have a significant risk of airflow obstruction. Am J Respir Crit Care Med. 2001;164:1419–1424. doi:10.1164/ajrccm.164.8.2105002

6. Hunninghake GM, Cho MH, Tesfaigzi Y, et al. MMP12, lung function, and COPD in high-risk populations. N Engl J Med. 2009;361(27):2599–2608. doi:10.1056/NEJMoa0904006

7. Ding Z, Wang K, Li J, et al. Association between glutathione S-transferase gene M1 and T1 polymorphisms and chronic obstructive pulmonary disease risk: a meta-analysis. Clin Genet. 2019;95(1):53–62. doi:10.1111/cge.13373

8. Cho MH, Boutaoui N, Klanderman BJ, et al. Variants in FAM13A are associated with chronic obstructive pulmonary disease. Nat Genet. 2010;42(3):200–202. doi:10.1038/ng.535

9. Pillai SG, Ge D, Zhu G, et al. A genome-wide association study in chronic obstructive pulmonary disease (COPD): identification of two major susceptibility loci. PLoS Genet. 2009;5(3):e1000421. doi:10.1371/journal.pgen.1000421

10. Soler Artigas M, Wain LV, Repapi E, et al. Effect of five genetic variants associated with lung function on the risk of chronic obstructive lung disease, and their joint effects on lung function. Am J Respir Crit Care Med. 2011;184(7):786–795. doi:10.1164/rccm.201102-0192OC

11. Repapi E, Sayers I, Wain LV, et al. Genome-wide association study identifies five loci associated with lung function. Nat Genet. 2010;42(1):36–44. doi:10.1038/ng.501

12. Han J, Lee Y, Yeom KH, et al. Molecular basis for the recognition of primary microRNAs by the Drosha-DGCR8 complex. Cell. 2006;125(5):887–901. doi:10.1016/j.cell.2006.03.043

13. Lund E, Guttinger S, Calado A, et al. Nuclear export of microRNA precursors. Science. 2004;303(15):95–98. doi:10.1126/science.1090599

14. Ambros V. The functions of animal microRNA s. Nature. 2004;431:350‑355. doi:10.1038/nature02871

15. Bartel DP. MicroRNAs: genomics, biogenesis, mechanism, and function. Cell. 2004;116:281‑297. doi:10.1016/s0092-8674(04)00045-5

16. Cao Z, Zhang N, Lou T, et al. microRNA-183 down-regulates the expression of BKCaβ1 protein that is related to the severity of chronic obstructive pulmonary disease. Hippokratia. 2014;18(4):328–332.

17. Donaldson A, Natanek SA, Lewis A, et al. Increased skeletal muscle-specific microRNA in the blood of patients with COPD. Thorax. 2013;68:1140–1149. doi:10.1136/thoraxjnl-2012-203129

18. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinform. 2008;9:559. doi:10.1186/1471-2105-9-559

19. Carlson MR, Zhang B, Fang Z, et al. Gene connectivity, function, and sequence conservation: predictions from modular yeast co-expression networks. BMC Genomics. 2006;7:40. doi:10.1186/1471-2164-7-40

20. Carter SL, Brechbuhler CM, Griffin M, et al. Gene co-expression network topology provides a framework for molecular characterization of cellular state. Bioinformatics. 2004;20(14):2242–2250. doi:10.1093/bioinformatics/bth234

21. Wan Q, Tang J, Han Y, et al. Co-expression modules construction by WGCNA and identify potential prognostic markers of uveal melanoma. Exp Eye Res. 2018;166:13–20. doi:10.1016/j.exer.2017.10.007

22. Zhai X, Xue Q, Liu Q, et al. Colon cancer recurrence-associated genes revealed by WGCNA coexpression network analysis. Mol Med Rep. 2017;16:6499–6505. doi:10.3892/mmr.2017.7412

23. Liu R, Zhang W, Liu ZQ, et al. Associating transcriptional modules with colon cancer survival through weighted gene co-expression network analysis. BMC Genomics. 2017;18:361. doi:10.1186/s12864-017-3761-z

24. Qin J, Yang T, Zeng N, et al. Differential coexpression networks in bronchiolitis and emphysema phenotypes reveal heterogeneous mechanisms of chronic obstructive pulmonary disease. J Cell Mol Med. 2019;23(10):6989–6999. doi:10.1111/jcmm.14585

25. Savarimuthu Francis SM, Davidson MR, Tan ME, et al. MicroRNA-34c is associated with emphysema severity and modulates SERPINE1 expression. BMC Genomics. 2014;15:88. doi:10.1186/1471-2164-15-88

26. Francis SM, Larsen JE, Pavey SJ, et al. Expression profiling identifies genes involved in emphysema severity. Respir Res. 2009;10(1):81. doi:10.1186/1465-9921-10-81

27. Langfelder P, Horvath S. Fast R functions for robust correlations and hierarchical clustering. J Stat Softw. 2012;46(11):11.

28. Xu P, Yang J, Liu J, et al. Identification of glioblastoma gene prognosis modules based on weighted gene co-expression network analysis. BMC Med Genomics. 2018;11(1):96. doi:10.1186/s12920-018-0407-1

29. Langfelder P, Horvath S. Eigengene networks for studying the relationships between co-expression modules. BMC Syst Biol. 2007;1:54. doi:10.1186/1752-0509-1-54

30. Zhang B, Horvath S. A general framework for weighted gene co-expression network analysis. Stat Appl Genet Mol Biol. 2005;4:17. doi:10.2202/1544-6115.1128

31. Russell SA, Bashaw GJ. Axon guidance pathways and the control of gene expression. Dev Dyn. 2018;247(4):571–580. doi:10.1002/dvdy.24609

32. Li HS, Chen JH, Wu W, et al. Vertebrate slit, a secreted ligand for the transmembrane protein roundabout, is a repellent for olfactory bulb axons. Cell. 1999;96(6):807–818. doi:10.1016/s0092-8674(00)80591-7

33. Tole S, Mukovozov IM, Huang YW, et al. The axonal repellent, Slit2, inhibits directional migration of circulating neutrophils. J Leukoc Biol. 2009;86(6):1403–1415. doi:10.1189/jlb.0609391

34. Lin YZ, Zhong XN, Chen X, Liang Y, Zhang H, Zhu DL. Roundabout signaling pathway involved in the pathogenesis of COPD by integrative bioinformatics analysis. Int J Chron Obstruct Pulmon Dis. 2019;14:2145–2162. doi:10.2147/COPD.S216050

35. He CW, Liao CP, Pan CL. Wnt signalling in the development of axon, dendrites and synapses. Open Biol. 2018;8(10):180116. doi:10.1098/rsob.180116

36. Qu J, Yue L, Gao J, Yao H. Perspectives on Wnt signal pathway in the pathogenesis and therapeutics of chronic obstructive pulmonary disease. J Pharmacol Exp Ther. 2019;369(3):473–480. doi:10.1124/jpet.118.256222

37. Baghy K, Tátrai P, Regős E, Kovalszky I. Proteoglycans in liver cancer. World J Gastroenterol. 2016;22(1):379–393. doi:10.3748/wjg.v22.i1.379

38. Frevert C, Wight TN. Matrix proteoglycans. In: Laurent GJ, editor. The Encyclopedia of Respiratory Medicine. London, UK: Elsevier; 2006:184–187.

39. Gill S, Wight TN, Frevert CW. Proteoglycans: key regulators of pulmonary inflammation and the innate immune response to lung infection. Anat Rec (Hoboken). 2010;293(6):968–981. doi:10.1002/ar.21094

40. Taylor KR, Gallo RL. Glycosaminoglycans and their proteoglycans: host-associated molecular patterns for initiation and modulation of inflammation. FASEB J. 2006;20:9–22. doi:10.1096/fj.05-4682rev

41. Coombe DR. Biological implications of glycosaminoglycan interactions with haemopoietic cytokines. Immunol Cell Biol. 2008;86:598–607. doi:10.1038/icb.2008.49

42. Wight TN. Arterial remodeling in vascular disease: a key role for hyaluronan and versican. Front Biosci. 2008;13:4933–4937. doi:10.2741/3052

43. Takahashi A, Majumdar A, Parameswaran H, et al. Proteoglycans maintain lung stability in an elastase-treated mouse model of emphysema. Am J Respir Cell Mol Biol. 2014;51(1):26–33. doi:10.1165/rcmb.2013-0179OC

44. Kim EK, Choi EJ. Compromised MAPK signaling in human diseases: an update. Arch Toxicol. 2015;89(6):867–882. doi:10.1007/s00204-015-1472-2

45. Sun Y, Liu WZ, Liu T, et al. Signaling pathway of MAPK/ERK in cell proliferation, differentiation, migration, senescence, and apoptosis. J Recept Signal Transduct Res. 2015;35(6):600–604. doi:10.3109/10799893.2015.1030412

46. Szatmary Z, Garabedian MJ, Vilcek J. Inhibition of glucocorticoid receptor-mediated transcriptional activation by p38 mitogen-activated protein (MAP) kinase. J Biol Chem. 2004;279(42):43708–43715. doi:10.1074/jbc.M406568200

47. Banerjee A, Koziol-White C, Panettieri R Jr. p38 MAPK inhibitors, IKK2 inhibitors, and TNFα inhibitors in COPD. Curr Opin Pharmacol. 2012;12(3):287–292. doi:10.1016/j.coph.2012.01.016

48. De Petris L, Brandén E, Herrmann R, et al. Diagnostic and prognostic role of plasma levels of two forms of cytokeratin 18 in patients with non-small-cell lung cancer. Eur J Cancer. 2011;47(1):131–137. doi:10.1016/j.ejca.2010.08.006

49. Zhang B, Wang J, Liu W, et al. Cytokeratin 18 knockdown decreases cell migration and increases chemosensitivity in non-small cell lung cancer. J Cancer Res Clin Oncol. 2016;142(12):2479–2487. doi:10.1007/s00432-016-2253-x

50. Yi H, Ku NO. Intermediate filaments of the lung. Histochem Cell Biol. 2013;140:65e69. doi:10.1007/s00418-013-1105-x

51. Nahm DH, Lee YE, Yim EJ. Identification of cytokeratin 18 as a bronchial epithelial autoantigen associated with nonallergic asthma. Am J Respir Crit Care Med. 2002;165:1536e1539. doi:10.1164/rccm.200201-009OC

52. Hacker S, Lambers C, Pollreisz A, et al. Increased soluble serum markers caspase-cleaved cytokeratin-18, histones, and ST2 indicate apoptotic turnover and chronic immune response in COPD. J Clin Lab Anal. 2009;23(6):372–379. doi:10.1002/jcla.20348

53. Lichtenauer M, Zimmermann M, Nickl S, et al. Transient hypoxia leads to increased serum levels of heat shock protein-27, −70 and caspase-cleaved cytokeratin 18. Clin Lab. 2014;60(2):323–328. doi:10.7754/clin.lab.2013.130303

54. Wang Z, Li W, Guo Q, et al. Insulin-like growth Factor-1 signaling in lung development, and inflammatory lung diseases. Biomed Res Int. 2018;2018:6057589. doi:10.1155/2018/6057589

55. Graff JW, Powers LS, Dickson AM, et al. Cigarette smoking decreases global microRNA expression in human alveolar macrophages. PLoS One. 2012;7(8):e44066. doi:10.1371/journal.pone.0044066

56. Wei P, Xie Y, Abel PW, et al. Transforming growth factor (TGF)-beta1-induced miR-133a inhibits myofibroblast differentiation and pulmonary fibrosis. Cell Death Dis. 2019;10:670. doi:10.1038/s41419-019-1873-x

57. Zhao H, Guo Y, Sun Y, et al. miR-181a/b-5p ameliorates inflammatory response in monocrotaline-induced pulmonary arterial hypertension by targeting endocan. J Cell Physiol. 2020;235:4422–4433. doi:10.1002/jcp.29318

58. Duan FG, Wang MF, Cao YB, et al. MicroRNA-421 confers paclitaxel resistance by binding to the KEAP1 30UTR and predicts poor survival in non-small cell lung cancer. Cell Death Dis. 2019;10:1–14. doi:10.1038/s41419-019-2031-1

59. Shen W, Liu J, Zhao G, et al. Repression of toll-like receptor-4 by microRNA-149-3p is associated with smoking-related COPD. Int J Chron Obstruct Pulmon Dis. 2017;12:705–715. doi:10.2147/COPD.S128031

60. Di T, Yang Y, Fu C, et al. Let-7 mediated airway remodelling in chronic obstructive pulmonary disease via the regulation of IL-6. Eur J Clin Invest. 2020;9:e13425. doi:10.1111/eci.13425

61. Liu S, Liu M, Dong L. The clinical value of lncRNA MALAT1 and its targets miR-125b, miR-133, miR-146a, and miR-203 for predicting disease progression in chronic obstructive pulmonary disease patients. J Clin Lab Anal. 2020;34(9):e23410. doi:10.1002/jcla.23410

62. Shigemura M, Lecuona E, Angulo M, et al. Hypercapnia increases airway smooth muscle contractility via caspase-7-mediated miR-133a-RhoA signaling. Sci Transl Med. 2018;10(457):1662. doi:10.1126/scitranslmed.aat1662

63. Pihtili A, Bingol Z, Kiyan E. Serum endocan levels in patients with stable COPD. Int J Chron Obstruct Pulmon Dis. 2018;13:3367–3372. doi:10.2147/COPD.S182731

64. In E, Kuluöztürk M, Turgut T, et al. Endocan as a potential biomarker of disease severity and exacerbations in COPD. Clin Respir J. 2021;15(4):445–453. doi:10.1111/crj.13320

65. Dai L, He J, Chen J, et al. The association of elevated circulating endocan levels with lung function decline in COPD patients. Int J Chron Obstruct Pulmon Dis. 2018;13:3699–3706. doi:10.2147/COPD.S175461

66. Yuan HS, Xiong DQ, Huang F, Cui J, Luo H. MicroRNA-421 inhibition alleviates bronchopulmonary dysplasia in a mouse model via targeting Fgf10. J Cell Biochem. 2019;120:16876–16887. doi:10.1002/jcb.28945

67. Barnes PJ. Oxidative stress-based therapeutics in COPD. Redox Biol. 2020;33:101544. doi:10.1016/j.redox.2020.101544

68. Jing X, Luan Z, Liu B. miR-558 reduces the damage of HBE cells exposed to cigarette smoke extract by targeting TNFRSF1A and inactivating TAK1/MAPK/NF-κB pathway. Immunol Invest. 2021;1–15. doi:10.1080/08820139.2021.1874977

69. Li X, Feng Y, Yang B, et al. A novel circular RNA, hsa_circ_0030998 suppresses lung cancer tumorigenesis and Taxol resistance by sponging miR-558. Mol Oncol. 2021;15(8):2235–2248. doi:10.1002/1878-0261.12852

70. Mao Y, He JX, Zhu M, et al. Circ0001320 inhibits lung cancer cell growth and invasion by regulating TNFAIP1 and TPM1 expression through sponging miR-558. Hum Cell. 2021;34(2):468–477. doi:10.1007/s13577-020-00453-4

71. Li XJ, Chen LW, Gao P, et al. MiR-587 acts as an oncogene in non-small-cell lung carcinoma via reducing CYLD expression. Eur Rev Med Pharmacol Sci. 2020;24(24):12741–12747. doi:10.26355/eurrev_202012_24173

72. Tian H, Wang X, Lu J, et al. MicroRNA-621 inhibits cell proliferation and metastasis in bladder cancer by suppressing Wnt/β-catenin signaling. Chem Biol Interact. 2019;308:244–251. doi:10.1016/j.cbi.2019.05.042

73. Wang YH, Yin YW, Zhou H, et al. miR-639 is associated with advanced cancer stages and promotes proliferation and migration of nasopharyngeal carcinoma. Oncol Lett. 2018;16(6):6903–6909. doi:10.3892/ol.2018.9512

74. Zhang Y, Zhang S, Yin J, Xu R. MiR-566 mediates cell migration and invasion in colon cancer cells by direct targeting of PSKH1. Cancer Cell Int. 2019;19:333. doi:10.1186/s12935-019-1053-1

75. Dai W, Liu S, Zhang J, et al. Vorinostat triggers miR-769-5p/3p-mediated suppression of proliferation and induces apoptosis via the STAT3-IGF1R-HDAC3 complex in human gastric cancer. Cancer Lett. 2021;521:196–209. doi:10.1016/j.canlet.2021.09.001

76. Xie L, Jiang T, Cheng A, et al. MiR-597 targeting 14-3-3σ enhances cellular invasion and EMT in nasopharyngeal carcinoma cells. Curr Mol Pharmacol. 2019;12(2):105–114. doi:10.2174/1874467212666181218113930

Creative Commons License © 2021 The Author(s). This work is published and licensed by Dove Medical Press Limited. The full terms of this license are available at and incorporate the Creative Commons Attribution - Non Commercial (unported, v3.0) License. By accessing the work you hereby accept the Terms. Non-commercial uses of the work are permitted without any further permission from Dove Medical Press Limited, provided the work is properly attributed. For permission for commercial use of this work, please see paragraphs 4.2 and 5 of our Terms.