Bioinformatics analyses of the pathogenesis and new biomarkers of chronic obstructive pulmonary disease.

BACKGROUND: Chronic obstructive pulmonary disease (COPD) is one of the major cause of global death. The purpose of our analysis was to detect a more reliable biomarker and small-molecule drug candidates and to identify the precise mechanisms involved in COPD.
METHODS: Three data sets were downloaded from the Gene Expression Omnibus database and analysed by Gene Expression Omnibus 2R. Functional enrichment analyses were performed by Metascape. We use the STRING data to build a protein-protein interaction network. The targets of differentially expressed microRNA (DE miRNA) were predicted by the miRWalk database. Small-molecule drugs were predicted on connectivity map.
RESULTS: A total of 181 differentially expressed genes and 35 DE miRNAs were confirmed. The protein-protein interaction network including all integrated differentially expressed genes was constructed, and 4 modules were filtrated. The module genes were relative to immune, inflammatory and oxidative stress functions according to a pathway analysis. The top 20 key genes were screened. Among the DE miRNAs found to be regulating key genes, miR-194-3p, MiR-502-5p, MiR-5088-5p, MiR-3127-5p, and miR-23a-5p might be the most significant due to their high number of connecting nodes in COPD. In addition, cephaeline, emetine, gabapentin, and amrinone were found to be potential drugs to treat COPD patients.
CONCLUSION: Our study suggests that miR-194-3p, miR-502-5p, and miR-23a-5p might participate in the nosogenesis of COPD. In addition, 4 potential small-molecule drugs were considered potentially useful for treating COPD patients.

Introduction

Chronic obstructive pulmonary disease (COPD) is a common, treatable chronic respiratory disease for the characteristics of continuous airflow restriction and repeated respiratory symptoms, which can developed into serious disease such as pulmonary heart disease and/or respiratory failure.[ As a common and frequent chronic respiratory disease, cases of COPD is gradually increasing, reaching 13.1% in 2015.[ At the same time, COPD is the eighth leading cause of disability worldwide,[ causing a huge economic burden to individuals and society. Furthermore, according to the World Health Organization, COPD is the fourth major cause of global death and is foretelled to be the third in the near future.[ Therefore, it is urgent to solve the problem of COPD. Although the pathogenesis has been repeatedly explored by countless people, the exact molecular mechanisms of COPD have not yet been fully clarified. Currently, we know that cigarette smoking and exposure to harmful particles are closely associated with COPD,[ and many smokers continue to experience disease-related changes after smoking cessation.[ However, not all smokers develop COPD. The pathogenesis of COPD might involve multiple risk factors, such as environment and heredity.[ In addition, effective, preventive and therapeutic measures are inadequate in COPD patients. Thus, we urgently need to further explore the molecular pathogenesis of this disease to find new and effective treatments.

MicroRNAs (miRNAs) are a group of approximately 19 to 25 nucleotides in length, noncoding RNAs that play critical roles in the regulation of basic biological activities such as proliferation, migration and apoptosis.[ Notably, some miRNAs have been identified to participate in the occurrence and development of COPD. Recently, miRNA-378 was thought to prevent the inflammation and thus prevent COPD development.[ In addition, miRNA-21 could play crucial role in the process of airway remodelling in COPD and is expected to become a possible therapeutic target for the treatment of COPD.[

In recent years, microarray and high-throughput sequencing, as emerging technologies, have become popular, and these technologies have gradually been developed based on big data integration and bioinformatics, which are leveraged to screen the hub genes related to the nosogenesis and outcome of different diseases. We can accurately find differences in transcriptome levels between a disease group and control group by microarray analysis. Based on bioinformatics analysis, a gene expression microarray provides broad prospects for target prediction and development of targeted drugs at the molecular level.[ A large number of biomarkers have been found in cancer and other diseases by bioinformatics microarray analysis.[ However, differentially expressed genes (DEGs) that are inconsistent, or even contradictory, may be revealed in different studies, which might be the result of individual or regional differences or even different statistical methods. Therefore, to minimize the constraints associated with a single data set, we integrated microarray data from multiple data sets to identify molecules more reliably and explore more precise mechanisms to use in treating COPD patients.

Materials & methods

Microarray data

Three data sets (GSE60399, GSE94916, and GSE102915) were retrieved and downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) with the keyword “chronic obstructive pulmonary disease”. These 3 datasets were based on the sequencing of RNA extracted from human peripheral blood mononuclear cells. GSE60399 and GSE94916 were related to gene expression profiling. The GSE94916 data set included 6 COPD patients and 6 healthy people based on the GPL18451 platform. The platform of the GSE60399 data set, with 160 COPD patients and 60 healthy peoples, was GPL20844.[ The GSE102915 data set comprised miRNA expression profiles from 6 COPD patients and 6 healthy people, and its platform was GPL20712.[

Identification of DEGs and DE miRNAs

GEO2R (https://www.ncbi.nlm.nih.gov/geo/geo2r) is an online tool for acquiring GEO array data,[ which are used to analyse gene and miRNA expression profiling. Gene IDs were converted to corresponding standard gene names according to the corresponding platform data. Then, the duplicates were removed. The DEGs and differentially expressed (DE) miRNAs were selected on the basis of a threshold P-value <.05 and log2-fold change (FC) >1 or <−1. Volcano plots and heatmaps of these DEGs and DE miRNAs were generated by R.

Functional and pathway enrichment analysis

Metascape is an online portal providing interactome analysis and comprehensive analysis[ that are used to perform gene ontology (GO) analysis and pathway analysis for integrated DEGs with the Kyoto Encyclopedia of Genes and Genomes (KEGG) database. Biological processes, cellular components (CC), molecular functions, and KEGG pathways were confirmed by analysis with a threshold P-value <.05. The results were saved as a text file, and then, the ggplot2 and magrittr packages in R (https://cn.stringdb.org/) were used to create visuals for analyses of portions of the results.

Protein–protein interaction (PPI) network

Protein–protein interaction (PPI) networks were built in the STRING data (https://string-db.org/).[ The minimum required interaction score was 0.40. After removing the disconnected nodes in the network, the results were saved in a csv file. Cytoscape 3.7.0 software was used for further processing. The molecular complex detection plug-in was served as a tool to select gene modules with MCODE score ≥3.50. A new method, maximal clique centrality (MCC), was recently proposed to be better than the other methods; therefore, we selected hub genes with the cytoHubba plugin using MCC methods.[ Furthermore, the KEGG analysis of modules was done by Metascape.

Gene set enrichment analysis

In addition, gene set enrichment analysis (GSEA) 4.02 software was used to conduct a GSEA analysis with gene set c2 (cp.kegg.v.6.2.symbols.gmt). GSE94916, an expression profile based on an array involving 12 samples from the GEO database, was analysed as a data set. The number of permutations was “1000”, and the permutation type was “gene_type”. A P-value <.01, FDR <0.25 and NES <−2 or >2 were the criteria for statistically significance.

Prediction of DE miRNA targets and integration of DEGs

The upregulated and downregulated DEGs in the GSE60399 data set were intersected with the upregulated and downregulated DEGs in GSE94916 data set, respectively, to obtain integrated DEGs for further analysis. This work was performed by the Calculate and draw custom Venn diagrams tool. The targets for the DE miRNAs from the GSE61741 data set were predicted based on miRWalk (http://mirwalk.umm.uni-heidelberg.de/) results. The genes predicted in at least 2 databases were considered to be targets of the DE miRNAs. Then, Cytoscape software[ was used to build a miRNA–mRNA network.

Identification of small-molecule candidate drugs

The connectivity map is an online gene expression profile database.[ Human cells were treated with small-molecule drugs, and the DEGs after treatment were used to establish a database connected with these small-molecule drugs, gene expression levels and diseases. This data enabled us to quickly identify small-molecule drugs highly associated with disease and the possible mechanisms of action of the small-molecule drugs.

Ethical review

This study is a bioinformatics analysis based on the published GEO public dataset. Animal and human experiments are not involved. Therefore, approval by the ethics committee is not required.

Results

After screening by GEO2R with P-value <.05 and absolute log2 FC >1, a total of 1722 DEGs were confirmed in the GSE60399 data set, including 1249 upregulated DEGs and 473 downregulated DEGs, and 5873 DEGs were confirmed in the GSE94916 data set, which included 66 upregulated DEGs and 5807 downregulated DEGs. In addition, 35 DE miRNAs involving 25 upregulated and 10 downregulated DE miRNAs were selected under the same screening conditions from the GSE102915 data set. Volcano plots and heat maps of the DEGs and DE miRNAs are shown in Figure 1A. Among these DEGs, 181 genes (16 upregulated DEGs and 165 downregulated DEGs) were filtrated in both 2 databases (Fig. 1B), and the DEGs were integrated for in-depth analysis (Table 1).

Figure 1

Differentially expressed genes in 3 data sets and confirmation of integrated differentially expressed genes in terms of gene expression profiling. (A) Volcano plot and heatmap of differentially expressed RNAs in the GSE60399, GSE94916 and GSE102915 data set. In Volcano plot, the red, green and black dots represent genes that are upregulated, downregulated and not significantly differentially expressed, respectively. In heatmap, upregulated genes are shown in red, and downregulated genes are shown in blue; the rows represent discovered RNAs in dataset, the columns represent samples. (B) Venn diagram of the upregulated and downregulated genes in the GSE60399 and GSE94916 data sets. GSE60399 is shown in blue, GSE94916 is shown in red.

Table 1

Integrated differentially expressed genes in COPD.

DEGs	Gene names
Upregulated	SLC25A37, CEACAM6, ANXA3, DEFA4, TMEM45A, ORM1, MCEMP1, HP, HPR, MPO, ARG1, IGFBP2, CAMP, AQP1, BPI, CA1
Downregulated	LMAN1, LOC93622, MRPL1, BCL2, CD6, RHOH, ATP10A, NDRG2, HLA-DOA, TAF9, METTL5, THOC7, KLRG1, ITK, PTPN4, MYC, ZNF18, NELL2, FARS2, PDCD6, ANKH, TRIM4, CHD2, LAT, CD28, OBFC1, TMEM50B, SNAPC5, TRAF3IP3, CD8B, TBRG1, TCF20, LY9, USP39, PASK, TRIM5, ABLIM1, RPP38, KLRB1, ZAP70, PTPRCAP, IL12RB1, ID3, PNPLA4, RPS23, ZNF669, PGRMC2, PLEKHB1, BTN3A2, APBB1, NIPSNAP1, MAL, LEF1, TTC24, DDX26B, MRFAP1L1, MAN1C1, RUVBL1, MAD1L1, C1orf52, NMT2, MLLT3, CNOT10, ZNF184, ATG9B, IVD, UBASH3A, CLK1, GIMAP2, SIGIRR, TIMM9, GTF2F1, COPS8, SNRPN, CD48, TNRC6C, CD3G, IL7R, CD96, TIMM44, ZFR, TCF7, GMDS, ZNF579, KLHL3, PINX1, SCML4, CD40 LG, CYCS, AES, MTMR1, RIOK1, EXOSC2, POLR1E, LDHB, NUP88, PIGP, FXN, IL27RA, FAIM3, ASTE1, P2RY8, PEX26, CNOT7, ARHGAP12, PRKACB, THEM4, BCL11B, GPR18, TNFRSF25, CCR7, BIN1, BUB3, PPP1R13 L, GATA3, LDLRAP1, SYNCRIP, PIGY, BZW2, AK5, CUTC, CD22, LPIN1, CYFIP2, SATB1, LETMD1, CD1C, FBXL12, TOMM20, NOSIP, GPX7, FLT3 LG, RPA2, RPS4X, RPS27A, GZMK, CDK6, LCK, KLRK1, P2RX5, CD248, LTB, CD7, XRRA1, AGPAT4, PYHIN1, ZNF101, CD3D, TBP, DPH5, EBAG9, CD244, RORA, TBCC, ACADSB, CDK5RAP1, PLD4, PRKY, ITM2A, CD5, PORCN, RNF157, HKDC1, DENND2D, RPS6KB1

COPD = chronic obstructive pulmonary disease, DEGs = differentially expressed genes.

Functional enrichment analyses of the integrated DEGs

The enrichment of the integrated DEGs including GO terms and KEGG pathways was performed. The top 20 terms and 10 significant enrichment items in each respect are listed in Figure 2. In GO terms, the integrated DEGs were enriched in immunity, T cell receptor signaling pathway, regulation of cell cycle, signal transduction by p53 class mediator, inflammatory response, Wnt signaling pathway, apoptosis and oxidative stress in biological processes category; the integrated DEGs were enriched in side of membrane and other organisms according to CC category and in transcription factor binding, protein kinase activity and hydrolase activity according to the molecular functions category. In the KEGG pathway analysis, the integrated DEGs were enriched in the T cell receptor signaling pathway, cytokine–cytokine receptor interaction, NF-kappa B signaling pathway, PI3K-Akt signaling pathway, and Jak-STAT signaling pathway.

Figure 2

Functional enrichment of integrated differentially expressed genes. (A) Top 20 clusters with their representative enriched terms. (B) Results of the top 10 GO terms (BP, CC, and MF) and top 10 KEGG pathways. BP = biological process, CC = cellular component, KEGG = Kyoto Encyclopedia of Genes and Genomes, MF = molecular function.

PPI network and analysis of the modules

A PPI network involving integrated DEGs was built on STRING data. As shown in Figure 3A, 4 modules were selected from the network according to the molecular complex detection-based analysis of the entire PPI network. The KEGG pathway analysis showed that the module genes are majorly contacted with the T cell receptor signaling pathway, Th1 and Th2 cell differentiation, Th17 cell differentiation and the NF-kappa B signaling pathway (Fig. 3B). The top 20 hub genes selected by MCC were CD28, CD5, IL7R, ZAP70, LCK, ITK, CCR7, CD40LG, KLRK1, GATA3, CD3D, CD48, CD3G, CD7, CD1C, KLRB1, LAT, CD8B, LY9, and CD244.

Figure 3

(A) Protein–protein interaction network construction and module analysis. (B) Kyoto Encyclopedia of Genes and Genomes pathway analysis of the 4 modules.

To more comprehensively explore the potential functions of the DEGs, we performed KEGG analysis in GSEA 4.0.2 software. As shown in Figure 4, the GSEA results also showed that the DEGs were associated with the T cell receptor signaling pathway. Some key genes LCK, ZAP70, LAT, CD28, ITK, CD40LG, CD8B, CD3G, and CD3D are in this pathway.

Figure 4

Results of the gene set enrichment analysis. KEGG = Kyoto Encyclopedia of Genes and Genomes.

Construction of the miRNA–mRNA network

In addition, 15,380 targeted genes were predicted by 35 DE miRNAs in miRWalk, and 152 targets were presented as integrated DEGs. Thus, a miRNA–mRNA network involving 35 DE miRNAs and 152 DEGs was constructed using Cytoscape software (Fig. 5A). To foster intuitive understanding of the regulatory mechanism of the key genes, the miRNA–mRNA networks of the key genes were shown in Figure 5B. MiR-502-5p, miR-3151-3p, miR-3127-5p, miR-4476, miR-23a-5p, miR-500a-5p, miR5088-5p, miR-1287-5p, and miR-194-3p regulate key genes. Among these, MiR-502-5p, MiR-5088-5p, MiR-3127-5p, and miR-23a-5p have a higher number of connecting nodes; therefore, they might be more significant in COPD.

Figure 5

(A) miRNA–mRNA network. (B) The network of the miRNAs and key genes. Triangles represent differentially expressed miRNAs, and rectangles represent integrated differentially expressed genes. Red is upregulated, and green is downregulated.

Prediction of small-molecule drugs

According to the analysis, DEGs in 2 gene expression databases were consistently found for COPD patients and healthy people, and the connectivity map database was utilized to screen candidate small-molecule drugs. Ten small-molecule drugs were screened according to the P-value <.05 and absolute enrichment >0.80 criteria (Table 2). Among those selected, cephaeline, emetine, gabapentin, and amrinone were considered potential drugs treating COPD patients for their significant negative correlation.

Table 2

Results of the connectivity map analysis.

Rank	CMap name	Mean	N	Enrichment	P-value	Specificity	Percent non-null
1	Cephaeline	−0.608	5	−0.944	0	0.0121	100
2	Emetine	−0.646	4	−0.938	0	0.0118	100
3	Biotin	0.652	3	0.928	.00066	0	100
4	Gabapentin	−0.439	4	−0.84	.00117	0	100
5	Rotenone	0.54	4	0.829	.00127	0	100
6	Pseudopelletierine	0.543	4	0.812	.00237	0	100
7	Amrinone	−0.416	4	−0.812	.00243	0.0227	75
8	Quinostatin	0.643	2	0.956	.00338	0.0884	100
9	Imidurea	0.544	3	0.817	.01212	0.0394	100
10	Cromoglicic acid	0.522	2	0.86	.03972	0.0382	100

CMap = connectivity map.

Discussion

In our study, 2 mRNA and 1 miRNA expression profiles were analysed to find hub genes participated in the nosogenesis of COPD. Two original COPD mRNA expression profiles and 1 COPD miRNA expression data set were found in the GEO database, each of which included comparisons between peripheral blood mononuclear cell samples from COPD patients and normal volunteers. After using GEO2R and biological integration synthesis analysis with a P-value <.05 and absolute log2 FC >1 threshold, 181 DEGs and 35 DE miRNAs were filtered.

MiRNAs constitute a class of noncoding endogenous RNAs that bind to complementary targets on mRNA in the 3-untranslated region, thereby inhibiting or silencing the expression of the targeted mRNA by cutting it or repressing translation.[ Recently, some miRNAs were found to upregulate rather than inhibit the expression of their targeted mRNA,[ but these examples represent a minority. In our study, we identified targeted DE miRNA–mRNA using the miRWalk network. For the top 20 key genes, when targeted mRNAs are negatively regulated by miRNAs, IL7R, ITK, CD244, ZAP70, CD5, CD28, CD8B, and KLRK1 are the potential targets of miR-23a-5p. The targets of miR-3127-5p are IL7R, KLRK1, and CD40LG. The targets of miR-502-5p are ITK, IL7R, and CD1C. CD28, CD3G, and ZAP70 are the targets of miR-5088-5p. Among these miRNAs, all these are upregulated, but their targeted mRNAs are downregulated. In some cases, miRNAs can induce the expression of their targets. When miR-194-3p is downregulated, its targets, CD8B, CD40LG, CCR7, and LAT, are downregulated.

The miR-194 family could regulate cell migration, proliferation,[ and inflammation suppression.[ Recently, miR-194-3p was reported to inhibit the proliferation and migration of fibroblasts and human trophoblast cells.[ MiR-194-3p was suggested as a potential biological marker for COPD, but the mechanism of its action hasn’t been discussed.[ Chemokine receptor 7 (CCR7) has been proposed to have a key role in oxidative stress and inflammation.[ CCR7 is the target of miR-194-3p; therefore, we hypothesized that COPD might be caused by diminished miR-194-3p targeting of CCR7. CCR7 and miR-194-3p can be potential targets for anti-inflammation and antiabnormal oxidative stress in COPD patients.

In addition, the miR-23 family plays a key role in T cell immune responses.[ In our analysis, miR-23a-5p is upregulated in COPD patients. Upregulated miR-23a-5p promotes cell proliferation and migration.[ In addition, miR-23a-5p has a close correlation with the NF-kappa B signaling pathway in regulating inflammation and immunity.[ The hub genes CD8B, CD28, ITK, and ZAP70 are all enriched in the T cell receptor signaling pathway. Additionally, ZAP70, KLRK1, and CD244 are enriched in the natural killer cell-mediated cytotoxicity process. CD8B, CD28, CD244, ITK, KLRK1, and ZAP70 are the targets of miR-23a-5p. Thus, these findings support the contention that miR-23a-5p might be related to immunity. Although the molecular mechanism of COPD remains unclear, the corresponding immune response might offer a potential explanation. MiR-23a-5p may be a new target for use in the immunotherapy of COPD.

IL-7R, one of the targets of miR-502-5p, is not only enriched in the PI3K-AKT signaling pathway but also in the Jak-STAT signaling pathway. In our analysis, IL-7R is downregulated in COPD patients. IL-7R is confirmed to be downregulated both in the expression profile microarray of end-stage chronic respiratory diseases and in subsequent QPCR validation of the microarray.[ This finding is consistent with our other analyses. Although miR-502-5p has not been reported in COPD until now, miR-502-5p might participate in the mechanism of COPD by targeting IL-7R.

Finally, cephaeline, emetine, gabapentin, and amrinone are considered to be potential drugs for treating COPD patients, according to our analysis. Among these drugs, research on cephaeline and amrinone has had very limited in the last 5 years. Emetine is both a synergistic agent of anticancer drugs[ and an antiviral drug.[ Recently, emetine was found to be a potential drug for treating pulmonary hypertension.[ While emetine has not been studied in COPD, it might be a possible drug to treat COPD due to its anti-inflammatory effects in the lung. Gabapentin usually serves as a kind of adjuvant for treating postoperative pain because of its antihyper-algesic properties.[ However, it has recently been found to downregulate the expression of inflammatory cytokines,[ which supports our results that gabapentin might be a potential drug to treat COPD patients. However, whether these drugs can be used in the clinic needs to be confirmed by basic and clinical trials.

Of course, there are some limits to our analyses. First, they were all based on bioinformatics analysis of public databases and have not been verified by experiments. In addition, the sample size we analysed was limited. To fill this gap, we intend to verify our findings with in vitro and in vivo experiments.

Conclusions

In summary, we found that CCR7 and IL-7R were downregulated by downregulated miR-194-3p and upregulated miR-502-5p, respectively. In addition, miR-23a-5p might also be a potential biomarker for COPD in patients. In addition, cephaeline, emetine, gabapentin, and amrinone are considered potential drugs for COPD patients. Our study contributes to further study of the nosogenesis of COPD. At the same time, it provides new latent biomarkers for COPD and new small-molecule drugs.

Author contributions

Conceptualization: Xiqian Xing.

Data curation: Jie Liu, Xiao Li, Liqiong Zhang.

Formal analysis: Jie Liu, Xiao Li.

Funding acquisition: Jiao Yang, Xiqian Xing.

Investigation: Shuanglan Xu.

Project administration: Jie Liu.

Resources: Liqiong Zhang.

Software: Shuanglan Xu.

Supervision: Jihua Zhang.

Validation: Yi Zhang, Jiao Yang.

Visualization: Yi Zhang.

Writing – original draft: Jie Liu.

Writing – review & editing: Xiaochao Yu, Jiao Yang, Xiqian Xing.