Integrating molecular interactions and gene expression to identify biomarkers and network modules of chronic obstructive pulmonary disease

Abstract

BACKGROUND:

Chronic obstructive pulmonary disease (COPD) causes chronic obstructive conditions, chronic bronchitis, and emphysema, and is a major cause of death worldwide. Although several efforts for identifying biomarkers and pathways have been made, specific causal COPD mechanism remains unknown.

OBJECTIVE:

This study combined biological interaction data with gene expression data for a better understanding of the biological process and network module for COPD.

METHODS:

Using a sparse network-based method, we selected 49 genes from peripheral blood mononuclear cell expression data of 136 subjects, including 42 ex-smoking controls and 94 subjects with COPD.

RESULTS:

These 49 genes might influence biological processes and molecular functions related to COPD. For example, our result suggests that FoxO signaling may contribute to the atrophy of COPD peripheral muscle tissues via oxidative stress.

CONCLUSIONS:

Our approach enhances the existing understanding of COPD disease pathogenesis and predicts new genetic markers and pathways that may influence COPD pathogenesis.

1. Introduction

Chronic obstructive pulmonary disease (COPD) is the third leading cause of death globally. According to recent estimates, from 2010 to 2030, the number of COPD cases in developed countries will increase by over 150% [1]. COPD is a complex disease, so determining the genetic risk factors for the disease has been challenging. Several studies have been conducted to meet this challenge. Qiu et al. [2] examined the gene expression in COPD patients with related diseases by genome-wide association studies (GWAS) to determine the functional effects of known susceptible genes and find new genes associated with the disease. Sakornsakolpat et al. [3] identified 82 loci that may be associated with either COPD or population-based measures of lung function. However, the GWAS method often generates vast genome-wide “hits” and the cost to examine these “hits” is high. In another study, Bahr et al. [4] used multiple linear regression to adjust covariates to identify COPD candidate genes and pathways. They showed that differentially expressed selected pathways in COPD subjects included those that have a role in the inflammatory response of the immune system. According to Huang et al. [5], the secondary metabolic pathway, the xenobiotic metabolic pathway, and the cellular response to xenobiotic stimuli possibly contribute to the development and advancement of COPD.

Although various efforts have been made for identifying biomarkers and pathways, the specific causal COPD mechanism still remains unraveled. Most of the information generated in these studies is based on single-group data, wherein, the relationships among biological regulatory networks, protein-protein interactions, signaling pathways, and well-known genes in the statistical framework on which they are based have not been considered. In such a framework, the lack of biological information leads to the stability of prediction factors and reduces the predictive ability of the model [6]. Therefore, prior biological knowledge should be incorporated into the model to acquire more reliable and biologically significant results [7, 8, 9].

A few such studies include those of Li and Li [10], Chen et al. [11], and Wang et al. [12] who attempted to use biological knowledge. In their model, genome knowledge was encoded by a network, and its graph structure determined the specific relationship (edge) between genes (nodes). They then integrated the Laplace matrix into the penalty in linear, logistics, and Cox regression models. The network particularly represented many types of biological information such as the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways [7], the correlation between genes [9, 13], functional interaction network [14], or protein-protein interactions (PPIs) [15]. These models which utilize a priori biological knowledge are often called “network-based” methods.

The hypothesis that complex diseases such as COPD arise and develop due to interactions between several interrelated pathogenic genes, is supported by a growing body of evidence, indicating that the evaluation of the influence of any single variant is complicated [16]. In this study, we hypothesize that combining biological interaction information with gene expression data would provide better insight into biological processes and network modules for COPD. Therefore, we tried to unravel the biological process involved in COPD using a network-based method and the results were promising. This study enhanced our understanding of the network module of COPD and predicted the associated new candidate genes and pathways affecting its pathogenesis.

2. Method

We applied a Laplacian constraint approach [10] to integrate the biological network for the analysis of the gene expression data. Let a network $G=(V,E)$ , where $V$ represents the genes with $p$ dimension and $E$ represents the connections between genes. Assuming $u$ and $v$ are connected, then $e_{uv}=1$ , otherwise $e_{uv}=0$ . The standard Laplacian transform $L$ for $G$ is shown as:

$L_{uv}=\{\begin{array}{cc}1-w_{uv}/d_u\hfill & \text{if}u=v\text{and}{d}_u\ne 0,\hfill \\ -w_{uv}/\sqrt{d_u{d}_v}\hfill & \text{if}u\text{and}v\text{are adjacent}\hfill \\ 0\hfill & \text{otherwise}.\hfill \end{array}$

where $d_u$ and $d_v$ are the degrees of genes $u$ and $v$ , respectively. For any fixed non-negative $\lambda$ , the network-based model is:

$L(\lambda ,\beta )=l(\beta )+\lambda {\beta }^{T}L\beta .$(1)

As shown in Eq. (1), the first term represents the loss function. The second term is a network-based penalty, which forms as a standard Laplacian matrix that captures the biological interaction knowledge.

Equation (1) is ill-posed in high-dimensional applications when the number of genes $p$ is greater than the sample size $n$ . Then, regularization approaches are widely applied to address this issue of large $p$ and small $n$ [6, 14, 17, 18, 19, 20, 21, 22, 23]. When a regularization term is added to Eq. (1), the sparse network constrained regression can be represented as Eq. (2).

$L({\lambda }_1,{\lambda }_2,\beta )=l(\beta )+{\lambda }_{1}P(\beta )+{\lambda }_2{\beta }^{T}L\beta$(2)

where ${\lambda }_1$ and ${\lambda }_2$ are regularization variables that are responsible for balancing the tradeoffs between fit and complexity. A popular regularization technique, Lasso ( $L_1$ ) has the regularization term $P(\beta )={\sum }_{j=1}^p|{\beta }_j|$ . However, in the lasso method, ${\lambda }_1$ should be tuned very carefully because model $\beta$ bears substantial bias if it is too large, and model $\beta$ may not be sufficiently sparse if ${\lambda }_1$ is too small. To deal with this problem, the smoothly clipped absolute deviation (SCAD) [24] method was proposed, which can be expressed as:

$P_{\lambda ,\text{𝑆𝐶𝐴𝐷}}(\beta )=\{\begin{array}{cc}\lambda |\beta |,\hfill & \text{if}0⩽|\beta |<\lambda ,\hfill \\ -{\displaystyle \frac{{\beta }^2-2\alpha \lambda |\beta |+{\lambda }^2}{2(\alpha -1)}},\hfill & \text{if}\lambda ⩽|\beta |<\alpha \lambda ,\hfill \\ {\displaystyle \frac{(\alpha -1){\lambda }^2}{2}},\hfill & \text{otherwise}.\hfill \end{array}$(3)

The SCAD approach has several advantages over Lasso, such as better estimation accuracy and the oracle property. Consequently, we penalized the network-based method using the SCAD approach as we proposed earlier [14]. Finally, the method we used in this study was:

$L({\lambda }_1,{\lambda }_2,\beta )=l(\beta )+P_{{\lambda }_1,{\lambda }_2,\mathit{\text{SCAD-Net}}}(\beta )$(4)

Where

$P_{{\lambda }_1,{\lambda }_2,\mathit{\text{SCAD-Net}}}(\beta )=P_{{\lambda }_1,\text{𝑆𝐶𝐴𝐷}}(\beta )+{\lambda }_2{\displaystyle \sum _{1⩽i<k⩽p;}}{w}_{ik}{\left({\displaystyle \frac{{\beta }_i}{\sqrt{d_i}}}-{\displaystyle \frac{{\beta }_k}{\sqrt{d_k}}}\right)}^2,$(5)

and $l(\beta )$ is defined as a logistics regression model.

We used the coordinate descent algorithm to solve Eq. (4). For more detailed information, please refer to our earlier publication [14].

Algorithm:
Input: Training dataset $\{X_{n\times p},y_n\}$ , ${\lambda }_1$ , ${\lambda }_2$ and $L$ .
Output: Model parameter $\beta$
Step 1: Update $\beta {(i)}^t,i=1,\mathrm{\dots },p$ .
Step 2: Let $t\leftarrow t+1$ , if $t<E$ , then repeat Step 1.

3. Results

3.1 Data description

To identify the key biological process for COPD, we included 136 subjects with COPD (42 ex-smoking controls and 94 COPD subjects), in the study. Expression data from peripheral blood mononuclear cells from these subjects were collected [4].

We mapped the dataset to an official gene symbol, and calculated average expression levels for multiple probe-sets mapping to the same gene. BioGrid provides the biological interaction network $L$ , which includes 14,621 genes or proteins and 327,721 interactions. The gene expression data were integrated into the network $L$ , and 102,292 edges and 7024 genes remained.

3.2 Construction of the model and selection of biomarkers

We used tenfold cross-validation on multiple dimensions to find optimal regularization parameters of the model. With the estimated tuning parameters and all the training data, a classifier having 49 genes was constructed (Table 1) with a training classification error of 4.26%. Among all the cutoff points, the one with the highest sum of sensitivity and specificity was chosen.

Table 1

The selected 49 genes from chronic obstructive pulmonary disease gene expression data

Genes
MMP1	PAK2	SLC25A38	NSUN5
TGFB2	SOD2	BACH2	PMS2P1
CTC1	FAM13A	NUDT16L1	CYFIP2
FCMR	EFHC2	SNRPN	THAP7
GPCPD1	NUMA1	HMCES	TMEM134
FOXP1	VAPB	ITPKB	ZNF775
OFD1	PTEN	NMT2	TSTD1
USP20	BCL6	FBXO46	DNLZ
RPARP-AS1	PPARD	LONP1	UBIAD1
RAB43	MTERF4	PATZ1	IL21R
MMP12	AQP9	HIVEP2	RASGRP2
FIP1L1	NOSIP	SFMBT1	SDR39U1
IRAK3

Among these 49 genes, we observed some interesting findings. For example, smoke-induced emphysema has been implicated with matrix metalloproteinase (MMP)-1, a collagen degrader, and MMP-12, an elastin degrader, at least in animal models [25]. However, robust epidemiological data about serum MMPs in COPD are scarce, and there is a massive gap between experimental research and clinical epidemiology. MMPs genes may play a significant role in COPD, so with this perspective, we believe that our study could be a step towards bridge this gap. FoxP1 is a transcription factor important for the development of lung epithelial tissue. Recent data from the UK Biobank, ECLIPSE, and COPD Gene cohorts implicate genetic variants in the FOXP1 gene as important predictors of airflow limitation. Moreover, loss of FoxP1 protein increases endoplasmic reticulum stress in lung epithelial cells that may contribute to COPD development [26]. While the biological function of the FAM13A gene product is poorly understood. genetic variants in the FAM13A gene might determine susceptibility to COPD and lung cancer [27]. Iwona et al. [28] confirmed that the FAM13A variants are predisposed to increased susceptibility to COPD.

These findings imply that the selected genes may contribute to or act as a marker for the pathophysiology of COPD.

3.3 Brief biological analysis

We then examined the 49 genes by GO and KEGG enrichment analyses (Figs 1 and 2). The results of GO analysis show that the selected 49 genes are involved in 17 significant pathways ( $p<$ 0.05), including regulation of B cell apoptosis, negative regulation of the production of cytokines involved in immune response, response to metal ions, muscle cell proliferation, aging, collagen metabolism, mitotic spindle organization, regulation of defense response, muscle structure development, extrinsic apoptotic signaling pathway, perinuclear region of cytoplasm, Golgi membrane, mitochondrial transport, anion transmembrane transport, chromatin binding, cellular response to abiotic stimulus and nucleocytoplasmic transport. KEGG analysis suggested that the selected genes are enriched in FoxO signaling pathway, cancer-associated pathways, microRNAs in cancer, and MAPK signaling pathways.

The enriched pathways may have possible roles in COPD. For example, the combination of chronic bronchiolitis and emphysema leading to COPD causes chronic airflow limitation and their ratio varies from patient to patient. B cell-related genes have been identified as being more prevalent in COPD with emphysema than in bronchiolitis, according to Faner et al. [29]. Likewise, Tang et al. [30] found an increased number of $T_C$ 2-like cytokine-expressing cells in the lungs of COPD patients. The genes associated with the negative regulation of cytokine production involved in the immune response pathway might explain why lung eosinophilia occurs during a COPD exacerbation.

This information can help us understand COPD pathobiology better, leading to new therapeutic possibilities for COPD.

We conducted a PPI enrichment analysis of the selected genes and also performed the Molecular Complex Detection (MCODE) algorithm to identify densely connected network components from the constructed PPI network. The MCODE network result is presented in Fig. 3.

We then performed pathway analysis of the MCODE network, and the result is shown in Table 2.

Table 2

Pathway analysis of the molecular complex detection network

GOID	Description	Log10 ( $p$ )
hsa04068	FoxO signaling pathway	$-$ 4.6
hsa05200	Pathways in cancer	$-$ 3.9
hsa05206	MicroRNAs in cancer	$-$ 3.3

Figure 1.

The gene ontology (GO) enrichment analysis.

Figure 2.

The Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses.

Figure 3.

The molecular complex detection algorithm resulted in the construction of protein-protein interaction (PPI) network from the selected 49 genes.

Atrophic muscle from COPD patients is shown to have increased expression of proteolysis pathway markers, such as transcription factors FoxO1 and FoxO3 [31]. Evidently, skeletal muscle atrophy involves increased expression of two E3 ubiquitin ligases, MuRF1 and atrogin-1, under the control of the FoxO family of transcription factors, particularly, FoxO1, activating protein degradation via the ubiquitin/proteasome system [32]. This information suggests that FoxO signaling may contribute to the atrophy of COPD peripheral muscle tissues via oxidative stress.

Data in Tables 1 and 2 and Figs 1–3 suggest that the selected 49 genes in this study might influence some biological processes and molecular functions related to COPD.

4. Discussion

Known for its chronic obstructive conditions, chronic bronchitis, and emphysema, COPD is the leading cause of death for the largest number of people worldwide. Although many efforts to identify biomarkers and pathways have been made, specific causal COPD mechanism remains unknown. This study combines biological interaction data with gene expression data for an improved understanding of the biological process and network module in COPD. We used a sparse network-based method and selected 49 genes from peripheral blood mononuclear cell expression data of 136 subjects, including 42 ex-smoking controls and 94 subjects with COPD. The results show that these 49 genes possibly influence biological processes and molecular functions related to COPD. For example, one of these genes coding for FoxO signaling may contribute to the atrophy of COPD peripheral muscle tissues via oxidative stress.

5. Conclusion

To conclude, our approach advances understanding of COPD disease pathogenesis and predicts new potential genetic markers and pathways that influence COPD pathogenesis. One limitation of this study was the lack of an in-depth verification of the selected genes and network modules.