Immune Characteristics Analysis and Transcriptional Regulation Prediction Based on Gene Signatures of Chronic Obstructive Pulmonary Disease.

PURPOSE: The variation in inflammation in chronic obstructive pulmonary disease (COPD) between individuals is genetically determined. This study aimed to identify gene signatures of COPD through bioinformatics analysis based on multiple gene sets and explore their immune characteristics and transcriptional regulation mechanisms.
METHODS: Data from four microarrays were downloaded from the Gene Expression Omnibus database to screen differentially expressed genes (DEGs) between COPD patients and controls. Weighted gene co-expression network analysis was applied to identify trait-related modules and then select key module-related DEGs. The optimized gene set of signatures was obtained using the least absolute shrinkage and selection operator (LASSO) regression analysis. The CIBERSORT algorithm and Pearson correlation test were used to analyze the relationship between gene signatures and immune cells. Finally, public databases were used to predict the transcription factors (TFs) and upstream miRNAs.
RESULTS: A total of 127 DEGs in COPD were identified from the combined dataset. By considering the intersection of DEGs and genes in two trait-related modules, 83 key module-related DEGs were identified, which were mainly enriched in interleukin-related pathways. Seven-gene signatures, including MTHFD2, KANK3, GFPT2, PHLDA1, HS3ST2, FGG, and RPS4Y1, were further selected using the LASSO algorithm. These gene signatures showed the predictive potential for COPD risks and were significantly correlated with 18 types of immune cells. Finally, nine miRNAs and three TFs were predicted to target MTHFD2, GFPT2, PHLDA1, and FGG.
CONCLUSION: We proposed the seven-gene-signature to predict COPD risk and explored its potential immune characteristics and regulatory mechanisms.

Highlights

1. Key module-related DEGs in COPD were found to involve in interleukin-related pathways.

2. An optimized gene set, including seven-gene signatures, was identified using LASSO analysis.

3. The seven-gene signature was significantly correlated with 18 types of immune cells in COPD.

4. The upstream miRNAs and TFs of gene signatures were predicted.

Introduction

Chronic obstructive pulmonary disease (COPD) is a complex heterogeneous disease characterized by persistent respiratory symptoms and airflow limitation because of respiratory and/or alveolar abnormalities.1,2 A forced expiratory volume in 1s (FEV1)/forced vital capacity (FVC) less than 0.70 after using bronchodilators is a necessary condition for the diagnosis of COPD.3 It is estimated to cause more than three million deaths worldwide annually, affecting about one in 10 adults aged over 45.4,5 Although smoking is the leading cause of COPD, the unsatisfactory prevalence of COPD in non-smokers suggests that chronic inflammation in patients may be caused by different cellular and pathophysiological changes in different genetic backgrounds.6 Some studies believe that genetic risk factors contribute to the development of airflow restriction, and the huge variation in lung inflammation between individuals is also genetically determined.7 Therefore, effective biomarkers are needed to predict risk and guide personalized treatment of patients with COPD.

Biomarkers are of great significance in the diagnosis, prognosis, and treatment of diseases. In recent years, many clinically relevant COPD biomarkers, including plasma fibrinogen, CRP, IL-6, IL-8, total bilirubin, SAA, SP-D, CCSP-16, MMP-8, and MMP-9, have been developed using the candidate gene approach.8–12 However, the influence of a single biomarker is often limited in predicting the disease. It has been reported that multiple biomarkers are more powerful than a single biomarker, and their combination with clinical variables can better predict disease progression.13 Based on shared data from high-throughput sequencing and microarray technology, it is possible to observe changes in the expression of numerous genes in different samples at one time.14 Moreover, bioinformatics methods that integrate multiple datasets can help identify reliable and repeatable genetic signatures associated with disease development.15 Therefore, data from four microarrays were used to analyze differentially expressed genes (DEGs) in this study. Then, a weighted gene co-expression network analysis (WGCNA) and the least absolute shrinkage and selection operator (LASSO) regression analysis were performed to obtain the optimized gene set from the trait-related module, which was also shown to have a great predictive ability among COPD patients.

Additionally, inflammation in the lungs of patients with COPD is associated with abnormal immune responses, and both the innate and adaptive immune systems are involved, linked by the activation of dendritic cells.16 Furthermore, an increase in the number of macrophages, neutrophils, and T and B lymphocytes has also been observed in the COPD process, and they are recruited in the immune cycle. They can secrete proinflammatory factors, such as cytokines, chemokines, growth factors, and lipid mediators.17,18 Considering the key roles of these immune cells, we attempted to identify the changes in the immune microenvironment defined by the gene signatures selected above and explain the molecular regulation mechanisms of upstream transcription factors (TFs) and miRNAs. A flowchart of this study is presented in .

Methods

Data Acquisition and Process

All data were downloaded from the Gene Expression Omnibus database (GEO, ). Using COPD, lung, and smokers as keywords, the relevant data of humans in the last five years were searched, and a total of four microarrays were obtained, including GSE37768, GSE106986, GSE103174, and GSE76925. Among them, GSE37768 included 18 lung tissue samples from COPD patients and nine healthy smokers and non-smokers, which were detected by GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array; GSE106986 included 14 lung tissue samples from COPD patients and five healthy non-smoker samples, which were detected by GPL13497 Agilent-026652 Whole Human Genome Microarray 4×44K v2; GSE103174 included 37 lung tissue samples from COPD patients and 10 healthy non-smoker samples which were detected by GPL13667 [HG-U219] Affymetrix Human Genome U219 Array; GSE76925 included 111 COPD samples and 40 smoker control samples which were detected by GPL10558 Illumina HumanHT-12 V4.0 expression beadchip. Among them, GSE76925 was used for the external validation. The detailed clinical information of the samples from these four datasets is summarized in .

Data Pre-Processing

The standardized probe expression value matrices of the above datasets were downloaded from GEO, and annotation files were downloaded from the corresponding platforms. Probes that did not match the gene symbol were removed, and when different probes were mapped to the same gene, the mean value of the probes was taken. Based on GSE37768, GSE106986, and GSE103174, the ComBat function in R sva (version 3.34.0, )19,20 was used to remove the batch effect, and a combined dataset was obtained for further analysis.

DEGs Analysis

To analyze the expression differences between COPD samples and controls in the combined dataset, the Limma package (version 3.10.3, )21 was used. The p-value obtained after the t-test was corrected by multiple tests using the Benjamini & Hochberg (BH) method, and the results were then evaluated in terms of difference multiples and significance. Adjusted p-value < 0.05, and |log fold change (FC)| > 0.585 were set as standards to select DEGs.

WGCNA

Genes with an absolute deviation of median expression in the top 2000 were selected for the WGCNA algorithm (version 1.61, ).22 WGCNA defined the correlation matrix of gene co-expression and the power exponential linkage function of the gene network, and the power value was captured when the square value of the correlation coefficient reached 0.9 for the first time. Then, a hierarchical clustering tree was constructed by calculating the dissimilitude degree of different nodes. The module was screened with a minimum number of genes in each gene network set to 30 and a pruning cut height of 0.3. Finally, the trait-related module was identified by analyzing the correlation between the module and disease phenotype with standards of p < 0.05 and |r| > 0.3. By considering the intersection of DEGs and genes in trait-related modules, key module-related DEGs were screened out using the R package VennDiagram (version 1.6.20, ).23

Enrichment Analysis

Based on the key module-related DEGs obtained above, over-representation analysis in the functional profiling module of gprofiler online tools ()24 was performed for gene ontology (GO) functions25 and pathways26 enrichment analyses. GO functions included biological process (BP), cellular component (CC), and molecular function (MF), whereas pathways were enriched in the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Reactome databases. BH adjusted p < 0.05 was considered statistically significant.

Screening and Validation of Gene Signatures

The key module-related DEGs were then incorporated into the LASSO regression model using the glmnet package (version 4.0–2, ),27 and the gene signatures were obtained with 20-fold cross-validation. Gene signatures were further verified on the combined datasets and GSE76925 by constructing a support vector machine (SVM) model using the e1071 SVM package (version 1.7–4, ).28 Receiver operator characteristic (ROC) curves were also created to verify the accuracy of the model.

Immune Cell Infiltration Analysis

The CIBERSORT deconvolution algorithm ()29 was used to estimate the infiltration abundance of 22 types of immune cells in each sample of the combined datasets, with the parameters set as perm = 100 and QN = F. Then, the Wilcoxon test was used to analyze the differences in the abundance of immune cell infiltration between the COPD and control groups, and a violin plot was created using the vioplot in R package (version 0.3.5, ).30 Finally, the Pearson correlation coefficient of gene signatures and immune cell infiltration abundance was calculated, and relation pairs (p < 0.05) were selected.

Predicting Upstream TFs and miRNA of Gene Signatures

The relationships between miRNA and gene signatures were predicted using miRWalk 3.0 ().31 The relationship pairs appearing in both miRDB ()32 and miRTarBase ()33 databases with scores over 0.9 were selected to construct a miRNA regulatory network. The TH-mRNA relations were predicted in the TRRUST database (),34 and the TF-mRNA and miRNA-mRNA relation pairs were integrated into a TF-mRNA-miRNA network.

Results

Screening of DEGs Between COPD and Healthy Controls

To remove the batch effect between GSE37768, GSE106986, and GSE103174, a principal component analysis was performed, and a combined dataset was obtained with thresholds of |logFC| > 0.585 and adjusted p < 0.05. A total of 127 DEGs between COPD and control samples in the combined dataset were screened out, as shown in Figure 1A. Detailed information on these DEGs is presented in . Among them, 104 DEGs were up-regulated, and 23 DEGs were down-regulated (Figure 1B).

Figure 1

Screening of DEGs between COPD and healthy controls. (A) The heatmap shows the DEGs between COPD and control samples from the combined dataset of GSE37768, GSE106986, and GSE103174, with thresholds of |logFC| > 0.585 and adjusted p < 0.05. (B) Volcano plot showing up-regulated (red) and down-regulated (blue) DEGs. The top five up-regulated and down-regulated DEGs ranked by |logFC| are labeled in red and blue, respectively.

WGCNA to Identify Trait-Related Modules and Key Genes

The WGCNA algorithm was then used to select trait-related modules and key genes. To satisfy the premise of a scale-free network distribution, the value of the power (weight parameter) of the adjacency matrix was explored. The selection range of network construction parameters was first set, then the scale-free distribution topology matrix was calculated, and the graph was created as shown in Figure 2A. When the square value of the correlation coefficient reaches 0.9 for the first time, the value of power is 6. A systematic cluster tree was obtained by calculating the coefficient of dissimilarity between genes. Then, the dynamic mixing cutting method was used to divide the modules, and finally, eight disease-related modules were obtained and labeled with different colors (Figure 2B). The genes involved in each module are listed in . Thereafter, the correlation between modules and disease traits was evaluated by calculating the correlation coefficient between the gene significance of each module and the diseased trait. The results (Figure 2C) suggested that module turquoise and module yellow had the highest correlation with disease traits, and these two modules were also identified as trait-related modules. By considering the intersection of genes in trait-related modules and DEGs obtained above, a total of 83 common genes were identified as key module-related DEGs (Figure 2D, ). Among them, 16 DEGs in module yellow were down-regulated, and 67 DEGs in module turquoise were up-regulated.

Figure 2

WGCNA algorithm was conducted to identify trait-related modules and genes. (A) Selection of power, weight parameter of the adjacency matrix. The x-axis represents the soft threshold (power), whereas the y-axis represents the square of the correlation coefficient between log2k and log2p(k) in the network. The red line represents the standard line, where the square value of the correlation coefficient reaches 0.9. (B) Module partition tree. The upper diagram represents the systematic cluster tree, and the different colors in the lower dynamic tree cut represent different gene modules. (C) Correlation between each module and disease traits. Statistical significance was set at p < 0.05. (D) The Venn plot shows the key module-related DEGs by taking the intersection of genes in trait modules and DEGs in COPD.

Functions and Pathways Enrichment Analyses

Based on the key module-related DEGs, enrichment analyses were performed to investigate their functions and pathways they may be involved in. As a result, 626 GO-BP, 13 GO-MF, 22 KEGG pathways, and 12 Reactome pathways were enriched. The top five results of each term ranking by p-value are shown in Figure 3A and B. It can be observed that these key module-related DEGs were mainly enriched in GO-BP of the response of cytokine and cellular response to cytokine stimulus; GO-MF of cytokine and chemokine activities; KEGG pathways of IL-17 and TNF signaling pathways; and Reactome pathways of IL-4, IL-13, and IL-10 signaling pathways. The genes involved in these key functions and pathways are listed in .

Figure 3

Functions and pathways enrichment analyses. (A) Top5 GO functions enriched in BP and MF. (B) Top5 KEGG pathways and Reactome pathways were enriched. The results of each term were sorted using p-values.

Identification of Gene Signatures

LASSO regression analysis was further performed to select the optimized gene set, and the parameters of the LASSO model are shown in Figure 4A. As a result, MTHFD2, KANK3, GFPT2, PHLDA1, HS3ST2, FGG, and RPS4Y1 were included in the optimized gene set when the lambda value was close to lambda.1se, and were thus identified as gene signatures. The expression of these gene signatures was then validated in the combined dataset (Figure 4B), and the results suggested that they were all significantly differentially expressed between COPD and healthy samples. ROC curves were also created to evaluate the efficacy of the seven-gene-signature in predicting COPD, as shown in Figure 4C and D. The results suggested that the area under the curve (AUC) created for the combined dataset was 0.795, whereas the AUCs of GSE37768, GSE106986, GSE103174, and GSE76925 were all over 0.700, thereby indicating a good predictive ability of gene signatures in COPD.

Figure 4

Identification of gene signatures using LASSO regression analysis. (A) The optimized gene set with seven DEGs was selected as the gene signature when the lambda value was close to lambda.1se. (B) Validation of gene signatures in the combined dataset. ****p < 0.0001. (C) ROC curve created for the combined dataset. (D) The ROC curves of GSE37768, GSE106986, GSE103174, and GSE76925.

Correlation Analysis of Gene Signatures and Immune Cell Infiltration Abundance

As described in the Methods section, the infiltration abundances of 22 types of immune cells were calculated for each sample. The differences in infiltration abundance between the COPD and control samples were analyzed, as shown in Figure 5A. The results suggested that fractions of plasma cells, CD8 T cells, resting NK cells, activated dendritic cells, and activated mast cells were significantly different between the two groups (p < 0.05). Thereafter, the correlation coefficients between gene signatures and immune cells were calculated, and relationships with p < 0.05 were considered statistically significant. As shown in Figure 5B and Table 1, FGG, GFPT2, HS3ST2, KANK3, MTHFD2, PHLDA1, and RPS4Y1 were significantly correlated with six, eight, eight, six, eleven, ten, and four types of immune cells. In these 53 gene-cell relationships, plasma cells were significantly related to MTHFD2, GFPT2, PHLDA1, HS3ST2, and RPS4Y1; CD8 T cells were significantly related to six gene signatures except PHLDA1; resting NK cells were significantly related to six gene signatures except RPS4Y1; activated dendritic cells were significantly related to FGG and PHLDA1; and activated mast cells were significantly related to five gene signatures except KANK3 and RPS4Y1.

Figure 5

Correlation analysis of gene signatures and immune cells. (A) The violin plot shows differences in immune cell fractions between COPD and control samples. Yellow and blue bars represent the cases and controls, respectively. *0.01

Table 1

Correlations Between Gene Signatures and Immune Cells

mRNA	Cell	p value	Corvalue	mRNA	Cell	p value	Corvalue
MTHFD2	B cells naive	0.038221549	−0.216460785	PHLDA1	NK cells activated	0.044817725	−0.209712723
MTHFD2	Plasma cells	0.009992671	0.267323475	PHLDA1	Monocytes	0.000754106	−0.345115438
MTHFD2	CD8 T cells	0.000941017	−0.339194798	PHLDA1	M0 Macrophages	0.037218299	0.217570938
MTHFD2	Naive CD4 T cells	0.006624417	−0.281183967	PHLDA1	Macrophages M1	0.005473208	−0.287396293
MTHFD2	Resting CD4 memory T cells	0.028518527	−0.228426768	PHLDA1	Activated Dendritic cells	0.034230963	0.221031971
MTHFD2	Resting NK cells	0.000560183	−0.352886978	PHLDA1	Resting Mast cells	4.38E-05	−0.412522682
MTHFD2	Monocytes	0.018706668	−0.244748154	PHLDA1	Activated Mast cells	8.93E-06	0.444735204
MTHFD2	M0 Macrophages	0.000685588	0.347627373	PHLDA1	Eosinophils	5.25E-05	0.408673492
MTHFD2	Resting Mast cells	5.98E-06	−0.452386635	HS3ST2	Plasma cells	0.002037279	0.317579771
MTHFD2	Activated Mast cells	0.000474456	0.357144707	HS3ST2	CD8 T cells	0.003800332	−0.298905194
MTHFD2	Eosinophils	0.000343708	0.365243811	HS3ST2	Naive CD4 T cells	0.005409791	−0.28777121
KANK3	CD8 T cells	0.000158799	0.383821809	HS3ST2	Resting NK cells	0.000447508	−0.358629682
KANK3	Naive CD4 T cells	0.014838555	0.253307477	HS3ST2	Monocytes	0.000978526	−0.338137567
KANK3	T cells follicular helper	0.012267997	−0.260140422	HS3ST2	M0 Macrophages	0.000142468	0.38634609
KANK3	Resting NK cells	4.37E-05	0.412591179	HS3ST2	Activated Mast cells	0.003092543	0.305210891
KANK3	M0 Macrophages	0.005401572	−0.287820088	HS3ST2	Eosinophils	0.015255688	0.252297256
KANK3	M2 Macrophages	0.011212825	−0.263312516	FGG	CD8 T cells	0.000532053	−0.354213995
GFPT2	Plasma cells	0.010632912	0.265168573	FGG	Resting NK cells	0.047402242	−0.207292668
GFPT2	CD8 T cells	0.004697775	−0.292272558	FGG	M0 Macrophages	0.013793222	0.255951605
GFPT2	Resting NK cells	0.001579882	−0.324868899	FGG	Activated Dendritic cells	8.64E-06	0.44537548
GFPT2	Monocytes	0.004339791	−0.29477023	FGG	Resting Mast cells	0.005283033	−0.28853245
GFPT2	Resting Mast cells	0.00020413	−0.377900936	FGG	Activated Mast cells	0.000234667	0.374564735
GFPT2	Activated Mast cells	4.99E-06	0.455739163	RPS4Y1	Plasma cells	0.004600289	0.292935454
GFPT2	Eosinophils	0.042492799	0.211990686	RPS4Y1	CD8 T cells	0.046803247	−0.20784364
GFPT2	Neutrophils	0.012392823	0.259781061	RPS4Y1	Naive CD4 T cells	0.038170325	−0.216516873
PHLDA1	Plasma cells	3.85E-05	0.415301636	RPS4Y1	M0 Macrophages	0.045106053	0.209437098
PHLDA1	Resting NK cells	5.60E-05	−0.407246483

Notes: Corvalues were p values adjusted by Benjamini & Hochberg method.

Correlation analysis of gene signatures and immune cells. (A) The violin plot shows differences in immune cell fractions between COPD and control samples. Yellow and blue bars represent the cases and controls, respectively. *0.01