Causal Path of COPD Progression-Associated Genes in Different Biological Samples

Abstract

Chronic obstructive pulmonary disease (COPD) is a progressive inflammatory disease with pulmonary and extra-pulmonary complications. Due to the disease’s systemic nature, many investigations investigated the genetic alterations in various biological samples. We aimed to infer causal genes in COPD’s pathogenesis in different biological samples using elastic-net logistic regression and the Structural Equation Model. Samples of small airway epithelial cells, bronchoalveolar lavage macrophages, lung tissue biopsy, sputum, and blood samples were selected (135, 70, 235, 143, and 226 samples, respectively). Elastic-net Logistic Regression analysis was implemented to identify the most important genes involved in COPD progression. Thirty-three candidate genes were identified as essential factors in the pathogenesis of COPD and regulation of lung function. Recognized candidate genes in small airway epithelial (SAE) cells have the highest area under the ROC curve (AUC = 97%, SD = 3.9%). Our analysis indicates that macrophages and epithelial cells are more influential in COPD progression at the transcriptome level.

Keywords:

• COPD
• candidate genes
• structure equation model
• elastic-net logistic regression

Introduction

Chronic obstructive pulmonary disease (COPD) is a progressive inflammatory disease characterized by airway obstruction that is among the first three causes of death worldwide [1,2]. Clinical presentations of COPD include emphysema, small airway obstructions, and chronic bronchitis.

Although the precise underlying processes of COPD pathogenesis are still unknown, the range of long-term exposure to many types of xenobiotics, such as smoking, occupational hazardous chemichals, and hereditary background, is stated [3], with no effective therapy or medicine [4,5]. Cigarette smoking is the most significant risk factor for the development of COPD [3,6]. Nonetheless, more than 20% of COPD patients are not smokers or ex-smokers [7,8]. Furthermore, additional research has shown that smokes emitted by solid fuels used for cooking or heating are not connected with COPD ethiology [9].

As a systemic multi-organ disease, multiple causes could alter the nature and consequences of COPD at the genomic level. Despite various transcriptomics studies on COPD, the critical role of hidden gene signatures is not clarified in different clinical samples [10–13]. It seems that a comprehensive analysis of heterogeneous biological samples is mandatory to find novel candidate genes in systemic pathogenesis and progression of COPD. The progression of COPD severity depends on the interaction of different cells present in the lung microenvironment. Various techniques have been utilized to interpret the disease complexity and pathogenesis elements. Many studies have been conducted based on top-down and/or bottom-up strategies to identify key players and regulators of the pathogenesis process, but nearly none investigated the level of direct/indirect impact of different cells on the staging and progression process of COPD. Specific statistical classification and prediction approaches can satisfy the need for high throughput high-dimensional transcriptome datasets [14]. Microarray analysis studies permit statistical and mathematical approaches to understand the association between thousands of genes in a disease’s specific stages. One of the primary challenges in microarray analysis is identifying genes, or groups of genes, that are differentially expressed in a disease or at different stages. More recently, machine-based learning algorithms have increasingly gained attention in bioinformatics and biology research [15,16]. In contrast, (regularization) based regression models (e.g. elastic-net logistic regression) have been widely used in microarray analysis [17].

We hypothesized that the different biological samples are not equally involved in the progression of COPD at the transcriptomics level. Also, we questioned if it would be possible to investigate the direct/indirect interaction among them. Therefore, this project was designed to identify novel important genes within different biological samples to provide novel therapeutic targets in COPD. Here, we fitted an elastic-net logistic regression model on transcriptome datasets to overcome overfitting and multi-collinearity issues, as common problems arise in the high throughput data analysis while representing a sparse and interpretable model in high dimensional datasets [12,18]. Following the identification of novel effective genes, structural equation modeling (SEM) was utilized to assess simultaneous direct and indirect effects of candidate genes and biological samples on COPD progression.

Methods

Datasets and data processing

A systematic search was performed in the PubMed database with the search strategy of COPD OR chronic obstructive pulmonary disease) AND (epitheli* OR macrophage OR basophil OR biopsy OR sputum OR BAL OR Bronchoalveolar lavage OR blood OR PBMC OR Lung tissue) AND (microarray OR transcriptomics OR RNA-seq). Then, omics data repository databases were searched to enrich the archived studies. After removing the duplicates, the studies with more than fifty samples and different stages of COPD with publicized raw data were included. Finally, studies were selected based on mostly matching transcriptomics platforms to minimize the cross-platform biases. The demographic and historical data of the selected datasets are provided in supplementary file 1.

Raw transcriptomics data from small airway epithelial cells (SAE) (GSE20257) [19], alveolar macrophages (GSE13896) [20], lung tissue (GSE47460) [21–23], sputum (GSE22148) [24], and blood samples (GSE54837) [25] were downloaded from GEO database. Initially, raw expression data were combined by the R package of "merging", as the combined dataset for each of biological samples included healthy controls and the patients in the GOLD stages of COPD [1–4]. After normalization of the datasets, the batch effect was removed and statistical comparison was performed by "Affy", "SVA", and "Limma" R packages, respectively [26,27]. The false discovery rate was corrected using the Benjamini-Hochberg correction method. The cutoff of the adjusted p-value (<0.0001) or absolute fold-change >2 was applied for the selection of differentially expressed genes. Since there were no healthy controls in sputum and blood datasets, the average normalized/standardized expression value of healthy controls from other datasets was used. Then Sensitivity analysis was performed to assess the robustness of the results obtained by differential expression analysis.

Elastic-net penalized logistic regression

Elastic-net logistic regression as a regularization method can be used for gene selection [28]. This statistical learning method was introduced as a compromise between ridge and lasso penalties. Elastic-net penalty combines the strengths of both ridge and lasso [29]. Based on the simulation and experimental studies on the genomic selection using regularized linear regression models, the elastic net regularization was found to outperform ridge and extensions of LASSO in high dimensional data [12, 30]. This penalty selects groups of correlated genes. Besides, it possesses optimized predictive performance compared with LASSO and ridge in the transcriptomics data [12]. Elastic-net logistic regression methods were performed by the "glmnet" R package [31]. Candidate genes were identified based on the below equation: $$\widehat{{\mathrm{\beta }}_{\text{Elastic net}}}=\mathrm{ }{\arg \min }_{\mathrm{\beta }}[-\sum _{\mathrm{i}=1}^{\mathrm{n}}\{{\mathrm{y}}_{\mathrm{i}}\mathrm{l}\mathrm{n}({\mathrm{\pi }}_{\mathrm{i}})+\left(1-{\mathrm{y}}_{\mathrm{i}}\right)\mathrm{l}\mathrm{n}(1-{\mathrm{\pi }}_{\mathrm{i}})\}+\mathrm{\alpha }{\text{ λ}}_1\sum _{\mathrm{j}=1}^{\mathrm{p}}\left|{\mathrm{\beta }}_{\mathrm{j}}\right|+\left(1-\alpha \right)\mathrm{ }{\mathrm{\lambda }}_2\sum _{\mathrm{j}=1}^{\mathrm{p}}{\beta }_j^2]$$

$\underbrace{\beta }$ is a vector of gene expression parameters, and lambda is tuning parameters for LASSO and ridge penalties, respectively. Beta parameters were estimated by coordinate descent as an optimization algorithm [31]. The best estimates for lambda parameters were derived by K-fold cross-validation. Alpha as a hyper-parameter, controls the distribution between the LASSO and ridge penalties. In this study, alpha was considered as a fixed value of 0.5.

Also, the statistical comparison was performed by the "pROC" R package between areas under the Roc curves (AUCs) for identified candidate genes in all biological samples [32]. The identification of novel candidate genes was performed by the elastic-net logistic regression method based on the trend of gene expression data over gold stages of COPD, and then, enrichment was performed using the ClueGO plugin of Cytoscape software.

Structural equation modeling

Structural equation modeling (SEM) includes causal modeling, analysis of covariance structures, and latent variable models. This modeling has many advantages as a generalization of multivariate multiple regression compared with conventional regression. SEM can assess multiple regression equations simultaneously, allowing identifying the strength and sign of direct and indirect effects for complex causal diagrams [33,34].

Path standardized coefficients (β) were calculated as the effect size that is free of scale and is robust against heterogeneity in the integration of genomics datasets. Goodness of fit (GOF) indices (e.g. The Root Mean Square Error of Approximation (RMSEA)<0.08, Standardized Root Mean Square Residual (SRMR)<0.08) were applied for assessing the fitness of the model. SEM was performed using the "lavaan" R package [35]. Numerical estimation of path coefficients (β) was derived by an iterative maximum likelihood algorithm. Multiple regression equations in SEM are represented below: $$\{\begin{array}{c}Y={\beta }_1{\mathit{gene}}_1+{\beta }_2{\mathit{gene}}_2+\dots +{\beta }_p{\mathit{gene}}_p \mathit{condidate} \mathit{genes}\mathit{\text{ in }}\mathit{SAE} \mathit{cells} \mathit{samples} \\ Y={\beta }_1{\mathit{gene}}_1+{\beta }_2{\mathit{gene}}_2+\dots +{\beta }_p{\mathit{gene}}_p \mathit{condidate} \mathit{genes}\mathit{\text{ in }}\mathit{alveolar} \mathit{macrophages} \mathit{samples} \\ Y={\beta }_1{\mathit{gene}}_1+{\beta }_2{\mathit{gene}}_2+\dots +{\beta }_p{\mathit{gene}}_p \mathit{condidate} \mathit{genes}\mathit{\text{ in }}\mathit{lung} \mathit{tissue} \mathit{samples}\\ Y={\beta }_1{\mathit{gene}}_1+{\beta }_2{\mathit{gene}}_2+\dots +{\beta }_p{\mathit{gene}}_p \mathit{condidate} \mathit{genes}\mathit{\text{ in }}\mathit{sputum} \mathit{samples}\\ Y={\beta }_1{\mathit{gene}}_1+{\beta }_2{\mathit{gene}}_2+\dots +{\beta }_p{\mathit{gene}}_p \mathit{condidate} \mathit{genes}\mathit{\text{ in }}\mathit{blood} \mathit{samples}\end{array}$$

Y in each equation is the number of stages in the biological sample. Beta coefficient is the strength and sign of candidate genes involved in the progression of COPD (Y).

Cross-validation, stability, and accuracy

The repeated k-fold cross-validation by bootstrapping is a good strategy to reduce the high variability of cross-validation [36]. Sensitivity analysis was performed to check the robustness of the results against changes in the training sets [37]. In the present study, the algorithm split the data set using repeated random 100 times sub-sampling in 5-fold cross-validation, permuting the sample labels every time. Cross-validated performances were summarized by observed sensitivities and specificities with standard deviation (SD). Furthermore, the area under the Receiver Operator Characteristic (ROC) curve (AUC) was used to calculate the precision of performance of the classifiers [38,39]. We used the bootstrap technique, which enables us to predict the fit of a model to a hypothetical testing set when an independent dataset is not available. The bootstrap approach allows us to use a computer to mimic the process of obtaining new datasets so that we can estimate the variability of our estimate without generating additional samples. Rather than repeatedly obtaining an independent dataset, we instead obtain distinct datasets by repeatedly sampling observations from the original data sets with replacement. To validate the selected genes with previous studies, literature mining was performed in PubMed. Interactive cluster heatmap was applied by the "heatmaply" R package [40].

Results

Differential analysis of genes expression data

Removing batch effects and normalizing data, according to the differential expression analysis of COPD vs. healthy samples, 918 probes from SAE, 1942 probes from lung tissue, 134 probes from blood, 1074 probes from alveolar macrophages, and 5768 probes from sputum samples were identified as the differentially expressed genes (adjusted p-value < 0.0001 or absolute fold-change >2) (Table 1).

Table 1. The summarized data indicating the primary, qualified, and differentially expressed probes in each biological sample.

Biological samples	GEO_ID	Platform	GOLD Stage (Number)	Samples	N. probes	N. differential probes
Lung tissue	GSE47460	GPL6480 GPL14550	Healthy=(N = 108) II=( N = 100) III=( N = 37) IV=( N = 83)	N = 328	N = 42,545	N = 1,942
Small Airway Epithelial cells (SAE)	GSE20257	GPL 570	Healthy=(N = 112) I=( N = 9) II=( N = 12) III=( N = 2)	N = 135	N = 54,675	N = 932
Blood samples	GSE54837	GPL570	I=( N = 88) II=( N = 70) III=( N = 55) IV=( N = 13)	N = 226	N = 54,675	N = 134
Macrophages	GSE13896	GPL6947	Healthy=( N = 54) I=( N = 8) II=( N = 3)	N = 65	N = 54,675	N = 1,074
Sputum	GSE22148	GPL570	II=( N = 71) III=( N = 59) IV=( N = 13)	N = 143	N = 54,675	N = 5,768

Gene selection and prediction

After adjustment of the effects of confounding variables (e.g. age and smoking status) by using Elastic-net penalized logistic regression, the total number of 33 genes was selected as associated factors with COPD progression with AUCs, sensitivities, and specificities in each biological sample (Table 2). According to statistical comparisons of AUCs of selected genes in different biological samples, genes identified in SAE cells and macrophages performed significantly better to predict the disease progression/stage. However, the AUC of the candidate genes in SAE samples was not significantly different (p-value = 0.478) compared to the AUC in macrophages to predict the disease stage (Figure 1).

Figure 1. Statistical comparison of AUCs among different biological samples based on the identified genes (* indicated significant difference at the level of 0.05, ** indicated significant difference at the level of 0.01, ***indicated significant difference at the level of 0.001).

Table 2. Probes and corresponding 33 candidate genes by elastic-net penalized logistic regression model for the association between the genes with COPD progression.

	Gene symbol	Probe ID	Up/down upregulated	Tissue	Epithelium	Blood	Macrophage	Sputum
	CCR4	A_23_P72989	Up	100%	–	–	–	–
	ITK	A_23_P354151	Up	84%	–	–	–	–
Tissue	RPUSD2	A_23_P309850	Down	82%	–	–	–	–
	RAB11B	A_23_P67748	Down	78%	–	–	–	–
	OXNAD1	A_24_P927189	Up	68%	–	–	–	–
	GPR171	A_23_P253317	Up	62%	–	–	–	–
	BTBD19	1557049_at	Up	–	100%	–	–	–
	THSD4	222835_at	Down	–	96%	–	–	–
	PPP4R4	233002_at	Down	–	95%	–	–	–
	NRG1	206343_s_at	Up	–	90%	–	–	–
	DNM3	209839_at	Up	–	89%	–	–	–
Epithelium	ITGA6	201656_at	Down	–	84%	–	–	–
	CD109	226545_at	Down	–	77%	–	–	–
	UHRF1	225655_at	Down	–	76%	–	–	–
	CST6	206595_at	Down	–	75%	–	–	–
	EPHB2	209589_s_at	Up	–	70%	–	–	–
	CDKN2A	207039_at	Up	–	70%	–	–	–
	KIAA1199	212942_s_at	Up	–	67%	–	–	–
	RGS20	210138_at	Down	–	63%	–	–	–
	SH3RF2	243582_at	Down	–	62%	–	–	–
	MTHFSD	244734_at	Up	–	–	100%	–	–
	CLEC7A	1554406_a_at	Up	–	–	89%	–	–
Blood	VCAN	211571_s_at	Up	–	–	66%	–	–
	PTPN4	236935_at	Down	–	–	–	100%	–
	CCDC37	242615_at	Up	–	–	–	77%	–
	GABARAPL1	208869_s_at	Up	–	–	–	67%	–
Macrophage	ADAMTSL1	1552808_at	Down	–	–	–	64%	–
	ATOH8	1558706_a_at	Down	–	–	–	64%	–
	SSBP1	202591_s_at	Up	–	–	–	62%	–
	SRPX	204955_at	Up	–	–	–	61%	–
	CHRFAM7A/CHRNA7	210123_s_at	Up	–	–	–	–	100%
Sputum	HSPA4	208814_at	Down	–	–	–	–	63%
	CADM1	232767_at	Up	–	–	–	–	62%
	AUC (SD)			0.92 (0.035)	0.97 (0.039)	0.73 (0.051)	0.96 (0.052)	0.82 (0.064)
	Sensitivity (SD)			0.86 (0.073)	0.89 (0.138)	0.46 (0.079)	0.88 (0.222)	0.62 (0.122)
	Specificity (SD)			0.80 (0.097)	0.93 (0.057)	0.82 (0.023)	0.95 (0.064)	0.82 (0.044)

Functional enrichment classified the novel genes into five groups, including "Regulation of CoA-transferase activity", "Vacuole organization", "dendritic spine organization", and "Cell adhesion molecules" (supplementary file 2). The expression level of candidate genes was measured and graphed among all biological samples (Lung Tissue, SAE, Blood, Macrophage, and Sputum) in each healthy and all COPD stages (Figure 2).

Figure 2. (A) Expression level of CCR4, ITK, RPUSD2, RAB11B, OXNAD1, and GPR171 in lung tissue. (B) Expression level of BTBD19, THSD4, PPP4R4, NRG1, DNM3, ITGA6, D109, UHRF1, CST6, EPHB2, CDKN2A, KIAA1199, RGS20, and SH3RF2 in SAE. (C) Expression level of MTHFSD, CLEC7A, and VCAN in Blood. (D) Expression level of PTPN4, CCDC37, GABARAPL1, ADAMTSL1, ATOH8, SSBP1, and SRPX in macrophages. (E) Expression level of CHRFAM7A/CHRNA7, HSPA4, and CADM1 in sputum.

For plotting co-expression patterns of selected genes among the patients, heatmap with agglomerative hierarchical clustering were plotted (Figure 3). Co-expression pattern of the selected genes resulted in four major clusters in COPD patients, including (OXNAD1, CCR4, ITK, and GPR171), (ADAMTSL1, THSD4, PPP4R4, ITGA6), (BTBD19, EPHB2, CHRFAM7A, SSBP1, GABARAPL1, ATOH8, PTPN4, MTHFSD, CCDC37, NRG1, CADM1, CLEC7A, VCAN), and (KIAA1199, DNM3, SRPX, CDKN2A, RPUSD2, RAB11B, HSPA4, RGS20, SH3RF2, CST6, CD109, UHRF1) based on the similarity of correlation coefficients/co-expression pattern (Figure 3). Of these 33 genes, 24 have previously been reported in the literature to be associated with lung diseases, including COPD or other lung disorders (Table 3). THSD4, PPP4R4, CDKN2A, CADM1, and NRG1, which has previously been detected in GWAS studies to determine single nucleotide polymorphisms (SNPs) in COPD and asthma, were among the mentioned 24 genes (https://www.ebi.ac.uk/gwas/home) [41–43]. However, we identified nine genes that have not been previously reported in COPD and other lung diseases, including RPUSD2, RAB11B, BTBD19, DNM3, SH3RF2, MTHFSD, ATOH8, SRPX, and HSPA4 (Table 3). These genes may represent novel potential biomarkers in the diagnosis and prognosis of COPD. The functional protein interaction network for the selected genes is illustrated in Figure 3, based on the STRING database (Figure 4).

Figure 3. Spearman’s rank correlation, co-expression, matrix between the selected genes in the COPD patients: heatmap for hierarchical clustering the 33 candidate genes based on their pattern of gene expression.

Figure 4. Co-expression network for the 33 selected genes from GeneMania. The size of the nodes indicates the co-expression degree between genes. The thickness of the edges indicates the weight of co-expression. The red-yellow color of the nodes/edges indicates the degree/weight, respectively.

Table 3. Confirmation of the association of 33 selected genes with COPD/or lung function by literature reviewing in PubMed databank.

Gene symbol	Probe ID	Number of studies	Associated diseases	References (PMIDs)
CCR4	A_23_P72989	4	idiopathic pulmonary fibrosis, lung cancer, lung metastasis in breast cancer	11590382, 16095529, 17168792, 23915095
ITK	A_23_P354151	7	idiopathic pulmonary fibrosis, sarcoidosis, allergic lung disease	16630934, 15323564, 12734350, 1646075, 26628680, 25512530, 24089408,
RPUSD2	A_23_P309850	0	–	–
RAB11B	A_23_P67748	0	–	–
OXNAD1	A_24_P927189	1	lung cancer	24040438
GPR171	A_23_P253317	1	lung cancer	26760963
BTBD19	1557049_at	0	–	–
THSD4	222835_at	8	COPD, airway diseases, asthma, pulmonary fibrosis	27564456, 24286382, 23932459, 23541324, 23409998, 22461431, 21965014, 20010834
PPP4R4	233002_at	1	COPD	28170284
NRG1	206343_s_at	23	lung cancer, COPD	31382039, 30988082, 30694715, 30568455, 30268483, 30069312, 29959202, and etc.
DNM3	209839_at	0	–	–
ITGA6	201656_at	4	Lung Fibrosis, idiopathic pulmonary fibrosis, lung cancer, prostate, head and neck cancers	31396340, 30936924, 28701000, 27143927
CD109	226545_at	5	Lung cancer, prostate and breast carcinoma, squamous epithelium, tumor cells	29113239, 28191885, 24667143, 17922683, 15116102
UHRF1	225655_at	8	Lung cancer, human ovarian cancer tissues, gastric and breast cancers, metastasis in hepatocellular carcinoma, renal cell carcinoma	30528265, 30008828, 29516630, 28849055, 27437769, 26695082, 21351083, 20517312
CST6	206595_at	1	lung cancer	24398667
EPHB2	209589_s_at	2	Allergic Rhinitis and Asthma, lung disease	28231727, 10037197
CDKN2A	207039_at	28	Lung cancer, COPD, head and neck cancer, malignant mesothelioma	30178167, 28487787, 27987577, and etc.
KIAA1199	212942_s_at	2	Non-small-cell lung cancer, cancer cells	30478628, 28901311
RGS20	210138_at	1	Lung cancer	29872324
SH3RF2	243582_at	0	–	–
MTHFSD	244734_at	0	–	–
CLEC7A	1554406_a_at	2	pulmonary fibrosis, lung inflammatory response	27852745, 27473664
VCAN	211571_s_at	7	Lung cancer, breast cancer, lung metastasis, asthma	27895126, 27581786, 27513329, 25044411, 21742797, 22392539, 23202429
PTPN4	236935_at	1	Non-small cell lung cancer (NSCLC)	26951513
CCDC37	242615_at	2		26200272, 22011669
GABARAPL1	208869_s_at	1	Non-small cell lung cancer (NSCLC)	26356813
ADAMTSL1	1552808_at	1	Lung cancer	29207642
ATOH8	1558706_a_at	0	–	–
SSBP1	202591_s_at	1	Lung cancer	28638454
SRPX	204955_at	0	–	–
CHRFAM7A/CHRNA7	210123_s_at	5	Non-small cell lung cancer (NSCLC), lung tumor, smoking-related lung cancers	25407004, 28283678, 28978081, 31096457, 30282908
HSPA4	208814_at	0	–	–
CADM1	232767_at	18	Lung cancer, lung tumor, lung epithelial cell apoptosis, lung fibroblasts, lung tumorigenesis	31069869, 23620770, 22429880, 22429880, and etc.

Causal pathway of selected candidate genes

Fitting the path diagram of selected genes (Figure 5), the genes in SAE, lung tissue, and sputum had more significant direct effects on COPD progression, respectively. In contrast, the identified genes in blood samples had less significant direct and indirect effects on COPD progression. Based on the magnitude of indirect standardized path coefficient, the novel genes in macrophages, lung tissue, SAE cells, and sputum affected COPD progression significantly indirectly compared with blood samples (Table 4). All goodness of fit indices indicated that the model has an acceptable fit (RMSEA = 0.059, p-value < 0.05; SRMR = 0.051).

Figure 5. The path diagram: A diagram based on the connections between selected genes in the studied biological samples.

Table 4. Standardized direct, indirect, and total effects of selected genes in studied biological samples on COPD progression.

Predictor → outcome	Direct effect (p-value)	Indirect effect (p-value)	Total effect (p-value)
Tissue → progression of COPD	0.79 (0.019)	0.24 (0.041)	1.03 (0.001)
Epithelium → progression of COPD	−0.82 (<0.001)	−0.22 (0.028)	−1.04 (<0.001)
Blood → progression of COPD	0.70 (0.038)	0.14 (0.048)	0.84 (0.032)
Macrophage → progression of COPD	−0.77 (0.001)	−0.29 (0.001)	−1.08 (<0.001)
Sputum → progression of COPD	0.76 (0.009)	0.18 (0.032)	0.94 (0.021)

Discussion

COPD patients can currently be classified based on clinical tests, including spirometry or St George’s Respiratory Questionnaire (SGRQ), but the need for further prognosis mandates identifying novel biomarkers that highly correlate with the disease progression. Nowadays, the application of high-throughput techniques as well as systems biology and machine learning approaches deeply alter the vision of researchers toward identification of novel diagnostic and prognostic biomarkers. Different studies have been conducted on the transcriptome of various cell lines involved in COPD’s pathogenicity. In this study, previously reported influencing cell types were studied to identify the most impacting genes in the development of COPD staging. Heterogeneity of samples make the algorithms vulnerable to identify higher number, more pure, and less noisy differentially expressed genes (DEGs), which in case narrows the discovery route of novel predicting genes.

The penalized logistic regression model and SEM are used to identify novel genes and to assess the connection network of the selected genes, respectively. The studies show that the elastic net often performs better than ridge and LASSO for the model selection consistency and prediction accuracy in microarray data [17]. The strengths of this study include the application of modern computational methods such as the statistical learning, SEM, validation of all of the results by literature review, and repeated cross-validation method (repeated 5-CV). However, some of the limitations of this study were the small sample size of disease stages, the lack of external validation dataset. On the other hand, applied novel methodologies in this study are robust against outliers and heterogeneity in gene expression data according to the bootstrap and cross-validation, standardized path coefficient in SEM and sensitivity analysis.

Analysis determined that identified genes in SAE cells, and macrophages had a higher accuracy rate than other biological samples in prediction of the disease progression. According to our results, twenty-four out of 33 genes were previously reported to associate with COPD, lung function (FVC, FEV₁ or the FEV₁/FVC ratio), and other lung related complications. Nine novel genes were identified with no background in COPD studies, including RPUSD2, RAB11B, BTBD19, DNM3, SH3RF2, MTHFSD, ATOH8, SRPX, and HSPA4. Amongst all, RPUSD2, RAB11B, BTBD19, DNM3, and MTHFSD were the most important novel genes in the analyses that can be nominated as novel important genes involved in the prognosis of COPD staging progression.

Contributing with the highest number of novel genes, epithelial cells appear to be most influential in pathogenesis of the COPD. The second rank in number of predicting genes belongs to macrophages. Macrophages have long been reported for their pillar role in inflammatory diseases and tissue remodeling processes. Enrichment analysis of novel candidate genes of epithelial cells demonstrated high contribution in cancer progression processes, initiation of inflammation and healing processes, and immune cell recruitment. On the other hand, some candidate genes impact hydrolyzing enzymes that leads to accumulation of mucus in bronchioles.

Thereupon, epithelial cells play a central direct role in COPD pathogenesis. However, macrophages are more indirectly in charge of inflammatory circumstances of lung microenvironment. Overall, it must be considered that macrophages have the upper hand compared to epithelial cells when it comes to COPD progression. Both M1 and M2 subclasses of macrophages are critical for healing processes. M1 cells remove damaged cells and M2 cells prepare an environment for proliferation and remodeling.

Macrophage candidate genes have high co-express with TWIST1 and PRRX1 transcription factors, which can be classified as oncogenes. These genes are mostly involved in genome replication and repair, along with autophagy and mitophagy processes. Comparing up and down-regulated candidate genes in macrophages reveal that downregulated genes mostly belong to antioxidant or proliferation initiation, while upregulated genes mostly belong to oxidant or tension responses.

Most candidate genes in tissue samples mainly contribute to chemokine signaling pathways, immune systems activation, cell-cell adhesion, and different G-protein related signaling cascades.

The co-expression network of candidate genes reveals that mostly tightly correlate with the immune system signaling mediators. We previously conducted research on COPD that highlighted the importance of epithelial cells in the progression of the disease. Also, 17 novel genes were introduced to be associated with the pathogenesis of COPD. PRKAR2B, GAD1, LINC00930, and SLITRK6 were the most important genes, which are consistent with the current study [12].

Conclusions

These novel genes may provide the basis for the development of therapeutics in COPD and its associated morbidities in the future. It is hoped that further studies on this issue would identify novel genes as biomarkers to help diagnosis and prognosis in COPD.

Authors’ contributions

SM participated in the design of the study and wrote the manuscript. HB and JS performed the data gathering. AE and MA performed the statistical analysis. AA and SAJ conceived of the study, and participated in its design and coordination, and helped to draft the manuscript. MS edited the manuscript. All authors read and approved the final manuscript.

Acknowledgments

This work was supported by Baqiyatallah University of Medical Sciences, Tehran, Iran.

Disclosure statement

The authors declare that they have no competing interests.

Funding

The author(s) reported there is no funding associated with the work featured in this article.

Data availability statement

Data sharing not applicable to this article as no datasets were generated or analyzed during the current study.