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

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 eu⁢v=1, otherwise eu⁢v=0. The standard Laplacian transform L for G is shown as:

where du and dv are the degrees of genes u and v, respectively. For any fixed non-negative λ, the network-based model is:

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).

where λ1 and λ2 are regularization variables that are responsible for balancing the tradeoffs between fit and complexity. A popular regularization technique, Lasso (L1) has the regularization term P⁢(β)=∑j=1p|βj|. However, in the lasso method, λ1 should be tuned very carefully because model β bears substantial bias if it is too large, and model β may not be sufficiently sparse if λ1 is too small. To deal with this problem, the smoothly clipped absolute deviation (SCAD) [24] method was proposed, which can be expressed as:

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:

Where

and l⁢(β) 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].

3.Results

3.3Brief 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 TC2-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.