Next Article in Journal
Tetracycline Analogs Inhibit Osteoclast Differentiation by Suppressing MMP-9-Mediated Histone H3 Cleavage
Next Article in Special Issue
Human Cysteine Cathepsins Degrade Immunoglobulin G In Vitro in a Predictable Manner
Previous Article in Journal
Electrical Stimulation of the Mesencephalic Locomotor Region Has No Impact on Blood–Brain Barrier Alterations after Cerebral Photothrombosis in Rats
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Consistent Biomarkers and Related Pathogenesis Underlying Asthma Revealed by Systems Biology Approach

1
College of Life Sciences, Chongqing Normal University, Chongqing 401331, China
2
School of Biological Information, Chongqing University of Posts and Telecommunications, Chongqing 400065, China
3
College of Laboratory Medicine, Chongqing Medical University, Chongqing 400046, China
*
Author to whom correspondence should be addressed.
These authors contributed equally to this work.
Int. J. Mol. Sci. 2019, 20(16), 4037; https://doi.org/10.3390/ijms20164037
Submission received: 21 July 2019 / Revised: 14 August 2019 / Accepted: 17 August 2019 / Published: 19 August 2019
(This article belongs to the Special Issue In Silico Analyses: Translating and Making Sense of Omics Data)

Abstract

:
Asthma is a common chronic airway disease worldwide. Due to its clinical and genetic heterogeneity, the cellular and molecular processes in asthma are highly complex and relatively unknown. To discover novel biomarkers and the molecular mechanisms underlying asthma, several studies have been conducted by focusing on gene expression patterns in epithelium through microarray analysis. However, few robust specific biomarkers were identified and some inconsistent results were observed. Therefore, it is imperative to conduct a robust analysis to solve these problems. Herein, an integrated gene expression analysis of ten independent, publicly available microarray data of bronchial epithelial cells from 348 asthmatic patients and 208 healthy controls was performed. As a result, 78 up- and 75 down-regulated genes were identified in bronchial epithelium of asthmatics. Comprehensive functional enrichment and pathway analysis revealed that response to chemical stimulus, extracellular region, pathways in cancer, and arachidonic acid metabolism were the four most significantly enriched terms. In the protein-protein interaction network, three main communities associated with cytoskeleton, response to lipid, and regulation of response to stimulus were established, and the most highly ranked 6 hub genes (up-regulated CD44, KRT6A, CEACAM5, SERPINB2, and down-regulated LTF and MUC5B) were identified and should be considered as new biomarkers. Pathway cross-talk analysis highlights that signaling pathways mediated by IL-4/13 and transcription factor HIF-1α and FOXA1 play crucial roles in the pathogenesis of asthma. Interestingly, three chemicals, polyphenol catechin, antibiotic lomefloxacin, and natural alkaloid boldine, were predicted and may be potential drugs for asthma treatment. Taken together, our findings shed new light on the common molecular pathogenesis mechanisms of asthma and provide theoretical support for further clinical therapeutic studies.

Graphical Abstract

1. Introduction

Asthma is a chronic inflammatory disease of the airways, caused by genetic or environmental and lifestyle factors, which is characterized by airway hyper-responsiveness (AHR), inflammation, and variable airflow obstruction [1]. Despite recent progress in the development of anti-asthmatic medication, asthma is still a major public health problem in the world. Its prevalence, morbidity, and mortality are still increasing [2]. According to the Global Asthma Report 2018, asthma affects over 339 million people, resulting in more than 1000 deaths every day (http://globalasthmareport.org). By 2025, an additional 100 million people may develop asthma [3]. However, this may be an underestimation due to underdiagnosis. A recent national survey in China found that the prevalence rate was ~4.2% in adults over 20, and over 45.7 million people suffered from asthma [4]. The total cost of asthma exceeded ~$81 billion per year in the USA and ~72€ billion in Europe [5]. Therefore, asthma not only affects human health seriously but also leads to serious socio-economic problems.
As a T-lymphocyte-controlled airway disease, asthma is caused by inflammation, overproduction of mucus, and airway remodeling, which results in hyper-reactivity and airway obstruction. The airway epithelial cells act as the first physical and immunological barrier and play an extremely important role in controlling inflammatory, immune, and regenerative responses to microorganism infections, allergens, and environmental pollutants that contribute to asthma pathogenesis [6]. Epithelial cells express pattern recognition receptors that detect the stimuli and release chemokines and cytokines, thereby bridging the innate and adaptive immune system cells [7]. Therefore, a better understanding of the underlying molecular mechanisms of the epithelium’s function in maintaining the integrity of the airways and its dysfunction in asthma is crucial for exploring how asthma is initiated and perpetuated and is also helpful in selecting new biomarkers for guiding therapeutic decision in asthma treatment.
Gene expression analysis is becoming more important in diagnostic fields allowing the identification of novel biomarkers relevant to diseases. The high throughput microarray is an efficient approach that facilitates gene expression analysis, new drug target identification, and novel gene function prediction. Woodruff et al. identified 22 differentially expressed genes (DEGs) (e.g., CLCA1, POSTN, and SERPINB1) in asthmatics when compared to healthy controls [8]. To identify gene expression patterns in children with asthma, Kicic et al. found 1612 genes (including 764 up-regulated and 848 down-regulated genes) were significantly differentially expressed in the epithelium of atopic asthmatics based on microarray analysis. The genes involved in the response to wounding were significantly reduced in epithelial cells, causing deficient immune responses and wound repair [9]. These findings strongly support that a fundamental alteration in the epithelium contributes to the initiation and progression of asthma. To link gene expression profiles to clinical asthma phenotype, Modena et al. performed a microarray analysis for 155 subjects with asthma and healthy controls and identified 1384 differentially expressed genes in bronchial epithelial cells. Further analysis revealed that some genes associated with type 2 inflammation, neuronal function, WNT signaling, and actin cytoskeleton could be considered as potential biomarkers for asthma [10]. A recent study using Affymetrix HT HG-U133+ PM GeneChips found that IL-13 response genes (POSTN, SERPINB2, CLCA1), mast cell mediator genes (CPA3 and TPSAB1), and cystatin genes (CST1, CST2 and CST4) were overexpressed in the epithelium of asthmatics [11], which was consistent with their previous study suggesting that activated T cells may be driving neutrophilic inflammation.
Although microarray experiments have generated long lists of genes with altered expression in asthmatics, some inconsistent findings were observed. Therefore, it is necessary to apply a systematic approach to combine different publicly available datasets to explore shared molecular mechanisms with asthma. The powerful meta-analytic technique has been well established to study the shared biological signatures between related diseases and pathophysiological conditions by merging multiple omics datasets [12,13].
In this study, we selected ten eligible microarray datasets to identify genes associated with dysfunctions of the bronchial epithelium in asthmatics. Differentially expressed genes were identified and functional annotations were performed for significant genes by using the microarray data integrating and analyzing techniques. Furthermore, approaches from systems biology, such as enrichment analysis and network analysis, were adopted to look for a better understanding of the molecular mechanisms underlying asthma.

2. Results and Discussion

2.1. Ten Independent Datasets Meeting the Criteria in This Study

In the study, 191 microarray datasets were obtained by keyword searching (as of 20th December 2018). After filtering based on dataset inclusion and exclusion criteria, 10 asthma-related microarray datasets from six disparate platforms (Affymetrix: HG-U95Av2, HG-U133A, HG-U133_Plus_2, HT_HG-U133_ Plus_PM, and Agilent: 4 × 44K G4112F and SurePrint G3 Human GE v3 8 × 60K) were selected. The detailed information of these datasets and sample descriptions is shown in Table 1 and Supplementary Material S1, respectively.
To compile expression data for data integrating and subsequent investigation, each dataset was preprocessed, normalized, and then integrated for further analysis following the steps illustrated in Figure 1.

2.2. Large Microarray Dataset Generated from Ten Selective Datasets

2.2.1. Sample Quality Control and Microarray Data Preprocessing

After the quality control, 1 (GSE470), 2 (GSE4302), 2 (GSE18965), 11 (GSE41861), 4 (GSE63142), 3 (GSE64913), 7 (GSE89809), and 4 (GSE104468) samples did not meet the cut-off criteria in array quality metrics assessment and were excluded for subsequent analysis; the QC reports are presented in Supplementary Material S2. Thus, a total of 556 samples were used for further analysis, containing 348 asthmatic patients and 208 healthy subjects. Through microarray data preprocessing, 10 gene expression matrices (genes in rows and samples in columns) were generated. For these matrices, all gene names were replaced by Entrez gene identifiers (IDs).

2.2.2. Data Integration and Batch Correction

Based on the results of the data preprocessing, probe annotations, and gene matching, a total of 8324 genes from 556 samples shared by six microarray platforms were selected (Supplementary Figure S1). After data integrating, batch effects were corrected using the ComBat algorithm. The clustering dendrogram showed that each dataset was clearly separated from the others before the batch correction (Figure 2A). After the batch correction, the samples from all datasets were well intermixed (Figure 2B), suggesting that the batch-adjusted data were suitable for further analysis.

2.3. Identification of DEGs Involved in Pathogenesis of Asthma

Based on feature selection, we identified 153 differentially expressed genes including 78 up-regulated and 75 down-regulated genes across the integrated dataset with the fold change > 1.2 and FDR < 0.05 (Table 2 and Supplementary Material S3). To identify genes most closely related to asthma, more stringent criteria (fold change > 1.5 and FDR < 0.05) were applied and only 8 up-regulated, 2 down-regulated genes were obtained. The up-regulated genes included CEACAM5 (carcinoembryonic antigen), CLCA1 (calcium-activated chloride channel regulator 1), POSTN (periostin), CPA3 (carboxypeptidase A3), SERPINB2 (plasminogen activator inhibitor-2), KRT6A (keratin 6A), CD44 (hyaluronate receptor), and MUC5AC (mucin 5AC), while LTF (lactotransferrin) and MUC5B (mucin 5B) were down-regulated. Gene expression patterns of all DEGs are shown in the volcano plot (Figure 3).

2.4. Functional Annotations and Enrichment Analysis of DEGs

2.4.1. Gene Ontology Analysis of DEGs

Based on DAVID analysis, a total of 12 GO functional terms were enriched with adjusted p < 0.01 and a minimum gene count > 6 (Table 3). Among these dysregulated GO terms, response to chemical stimulus (GO:0042221, Z-score = −0.34) was the most significantly enriched into biological process (BP) (Adj. p = 2.64 × 10−4), while membrane fraction (GO:0005624, Z-score = −0.89), insoluble fraction (GO:0005626, Z-score = −0.89), and cell fraction (GO:0000267, Z-score = −0.82) from cellular components (CC) were considered as the top three decreased terms (Figure 4). Interestingly, two GO terms relevant to asthma were also enriched: increased secretory granule (GO:0030141, Z-score = 0.63) and more active extracellular region (GO:0005576, Z-score = 0.42).

2.4.2. KEGG Pathway Analysis of DEGs

Based on KEGG pathway mapping, DEGs were significantly enriched in 12 pathways (FDR < 0.05) involved in cancer, arachidonic acid metabolism, linoleic acid metabolism, calcium signaling, aldosterone-regulated sodium reabsorption, and so on (Table 4). It is worth noting that 7 down- and 3 up-regulated genes involved in cancer were also enriched (FDR = 3.24 × 10−5, Z-score = −1.26), which is consistent with the previous report that asthma status was associated with decreased risk of aggressive bladder cancer [22]. Moreover, the pathway associated with bladder cancer was also enriched (FDR = 1.17 × 10−2) with a suppressed trend (Z-score = −0.58). Arachidonic acid metabolism was the second most significant pathway (FDR = 1.42 × 10−4). Its abnormal metabolism has been observed in asthma [23,24]. Dysregulated linoleic acid metabolism (FDR = 7.69 × 10−3) was the third enriched pathway. Woods et al. proposed that fatty acid levels were associated with the risk of asthma in young adults [25]. Black et al. argued that increasing dietary intake of n-6 linoleic acid had resulted in increased arachidonic acid and PGE2 production, with a consequent increase in the likelihood of asthma [26]. Our results also found the calcium signaling pathway was significantly suppressed in asthma (FDR = 1.17 × 10−2, Z-score = −1.34). An abnormal calcium signaling pathway has been linked with many diseases. Mahn et al. reported that the dysregulation of Ca2+ homeostasis was probably related to abnormal asthmatic phenotype [27]. Our result supports the fact that peroxisome proliferator-activated receptors (PPAR) are associated with obstructive lung disease [28]. Thus, PPAR may also be a potential target of asthma. Interestingly, a novel pathway named Hematopoietic cell lineage was significantly enriched in this study, suggesting some components in this pathway may be associated with asthma pathogenesis.

2.4.3. Potential Target Sites of Transcription Factors and Regulatory MicroRNAs

Using the GSEA online tool, 9 DNA-binding motifs with FDR < 0.05 were identified and are listed in Supplementary Material S4. Among them, a forked transcription factor FOXO4 and its downstream targeting genes RIT1, CLC, FKBP5, AGR2, CD36, RUNX2, NTRK2, GRK5, GATA2, ST6GAL1, SLC7A1, SLC18A2, and SERPINB10 were enriched. The previous studies showed that FOXOs were involved in multiple cellular processes, such as stress resistance, cell cycle, apoptosis, and metabolism [29]. In asthmatic patients, the increased expression of AGR2 and RUNX2 would lead to excessive secretion of mucin [30,31]. GATA2 positively regulated the expression of IL-33 receptor (ST2), which stimulated the production of several pro-inflammatory factors in mast cells [32]. Furthermore, the increased expression of IL-13 promoted the sialylation of MUC4β N-glycans by ST6GAL1 and inhibited the proliferation of damaged epithelial cells in asthmatic patients [33]. The damaged epithelial cells can produce fibroblast-promoting growth factors, which aggravated the airway remodeling response [34].
Similarly, 24 potential regulatory microRNAs were identified (FDR < 0.05) (Supplementary Material S5). A previous study showed that the expression level of microRNA-181a was higher in the acute asthma patients compared to the healthy controls at the beginning of asthma, then dropped to the control level when there were no new airway stimuli. Therefore, microRNA-181a is a potential pro-inflammatory factor in asthma [35]. Mohamed et al. reported that the overexpression of microRNA-26a would increase hypertrophy of human airway smooth muscle cells and promote airway remodeling by inhibition of GSK-3β [36].

2.4.4. Effect of Chromosomal Position on the Expression of DEGs

For many chromosomal loci, gene mutation and dysregulation of gene expression often result in many diseases including prostate cancer [37] and breast cancer [38]. Genome-wide linkage analysis of asthma and airway responsiveness showed that chromosome 12q24.31 contained a locus that was crucial in intermediate phenotype for asthma in a Hispanic population of Costa Rica [39]. Subsequently, Ferreira et al. found that chromosome 11q13.5 locus was significantly associated with the risk of allergic sensitization, in turn, increasing the risk of subsequent development of asthma in Australia [40]. In this study, gene position enrichment analysis of 153 DEGs revealed that 9 asthma-related genes were located on chromosome 3 (Figure 5). Among them, 5 up-regulated genes (CSTA, GATA2, CPA3, P2RY14, and AADAC) and 1 down-regulated gene (HEG1) were located in q21 region (FDR = 9.65 × 10−5), while the remaining 3 down-regulated genes (LTF, ITPR1, BHLHE40) were located out of q21 region on chromosome 3. Our results indicated that chr3q21 was a new locus associated with asthma.

2.5. Identification of Hub Genes Based on PPI Network Construction

To better understand the mechanism of asthma and identify the crucial hub genes among the DEGs, PPI network was constructed using network analysis, which enables the analysis of protein-protein interaction for multiple genes using NetworkAnalyst [41]. To categorize gene functions, genes in this network were divided into three groups: hub genes highly interacting with other neighbor genes (CD44, KRT6A, LTF, SERPINB2, MUC5B, and CEACAM5), bridge genes connecting different communities (e.g., MUC5AC, CLCA1, POSTN, SP1, and EGFR) and leaf-node genes with one neighbor (e.g., CPA3, BAG1, MNF1, and APP), as shown in Figure 6.
Within this network, CD44 encoding a cell surface adhesion receptor was the most highly ranked hub gene (degree = 59). A previous study showed that CD44 was significantly up-regulated in the bronchial epithelium in asthmatics compared with the controls [42]. Higher expression of CD44 in blood eosinophils could regulate the recruitment and function of leukocytes in asthmatics and are considered as a marker of bronchial asthma [43]. KRT6A are a member of type II epithelial keratins, which are intermediate filament-forming proteins that provide mechanical support and fulfil a variety of additional functions in epithelial cells. The presence of KRT6A was associated with pachyonychia congenital [44], as well as renal carcinoma [45] and breast cancer progression [46]. Genome-wide expression profiling with linkage analysis revealed that the transcription level of KRT6A significantly increased in peripheral blood mononuclear cells (PBMCs), airway brushing cells (ABCs), and bronchioalveolar lavage (BAL) when asthmatics were exposed to cockroach allergen [47]. In line with this, a recent proteomic study of asthma also found that KRT6A was up-regulated in the blood of asthmatics [48]. Therefore, KRT6A was considered as an important marker for asthma. Lactotransferrin (LTF), also called lactoferrin, is an iron-binding glycoprotein and servers as an immune-modulator and anti-inflammatory factor. Transcriptome analysis showed that LTF was significantly up-regulated during the development of asthma [49]. LTF could reduce pollen-induced airway inflammation in a mouse model of allergic asthma [50]. Immunoglobulin receptor CEACAM5, also known as CD66e, is a member of the carcinoembryonic antigen (CEA) family and involved in cell signaling, cell proliferation, cell repair responses, and the maintenance of the intact bronchial epithelium [51]. A previous study showed that the expression of CEACAM5 was increased in smoking and non-smoking severe asthma [52]. Consistent with this, the up-regulation of CEACAM5 was also observed in the epithelium in severe neutrophilic asthma [53].
Within bridge genes, MUC5AC is a glycoprotein and significantly increased when the respiratory tract was exposed to external stimulus [54]. High expression of MUC5AC was observed in bronchial epithelial cells of patients with severe asthma [55]. A recent study showed that IL-13 response genes (SERPINB2, POSTN, and CLCA1), protease genes (CPA3 and TPSAB1), and a nitric oxide synthase gene (NOS2) were significantly up-regulated in the epithelium of mild asthma [11].
To provide a comprehensive PPI network structure, the well-established network community recognition algorithm InfoMap [56] was applied and three distinct communities (marked with royal blue, dark turquoise or magenta areas) were obtained (p < 0.001). Combined with gene enrichment analysis, we found that they were associated with the regulation of response to stimulus, cytoskeleton, and response to lipid, respectively (Figure 6). This result suggested that the occurrence of asthma may be associated with the disturbance of lipid and abnormal cytoskeleton when epithelium cells were exposed to the stimulus.

2.6. Identification of Candidate Small Molecules

Using the 10 DEGs with fold change > 1.5 as seeds, 20 potential asthma-related chemical molecules with |score| > 0.6 and p < 0.05 were screened using a Connectivity Map (CMap)-based systems approach. Among them, 17 molecules probably promote the progression of asthma, such as Prestwick-1082, ricinine, and milrinone. On the contrary, 3 molecules have potential therapeutic effects, including catechin, lomefloxacin, and boldine (Table 5).
Catechins were one kind of polyphenols and showed various biological and pharmacological activities in antioxidative [57], anti-carcinogenic [58], and antiallergic effects [59]. Catechins from green tea could significantly inhibit the migration of inflammatory cells by suppressing MMP-9 expression and ROS generation in endothelial cells in a murine model of asthma [60]. Recently, Patel et al. confirmed that catechins from Acacia catechu showed inhibitory effects on ovalbumin-induced allergic asthma by inhibiting the activity of the histidine decarboxylase enzyme, and the author suggested that catechin can be used for the treatment of allergic inflammatory disease in humans [61].
Grassi et al. found that antibiotic lomefloxacin can treat acute exacerbations of chronic bronchitis [62]. Using the signalome screening method, Todd et al. reported that lomefloxacin showed significant effects on multiple pro-inflammatory signaling pathways that are constitutively activated in the auto-inflammatory disease TNF receptor [63]. However, no study has been performed concerning the usage of lomefloxacin for asthma treatment.
Boldine is a natural alkaloid from the leaves and bark of the Chilean boldo tree, Peumus boldus, and had anti-inflammatory [64] and antioxidant properties [65]. However, potential therapeutic effects on asthma have not been reported.
Taken together, the compounds catechin, lomefloxacin, and boldine may be potential drugs for the treatment of asthma caused by inflammation and allergy. However, further experiments will be needed to confirm the meta-analysis results.

2.7. Crosstalk Pathway of Asthma

To further explore the molecular mechanism of disease pathogenesis, functional enrichment analysis based on the BIOCARTA, REACTOME, and Pathway Interaction Database was performed by using the ToppGene Suite with the threshold of BH-FDR < 0.05.
Six categories were enriched, including 1) genes encoding ECM and ECM-associated proteins (M5889); 2) ECM-affiliated proteins, ECM-regulators and secreted factors (M5885); 3) HIF-1α transcription factor network; 4) FOXA1 transcription factor network; 5) termination of O-glycan biosynthesis, and 6) interleukin-4 and 13 signaling (Table 6). Interestingly, two asthma-associated pathways (M5889 and M5885) were enriched from the BIOCARTA database. Changes observed during airway remodeling in chronic asthmatic patients include excessive extracellular matrix (ECM) production and collagen deposition, increased airway smooth muscle mass, and mucus hypersecretion [66]. Abnormal deposition of ECM protein and/or ECM-associated proteins causes airway stiffening and narrowing, and differences in ECM protein expression may represent a specific asthma phenotype [67]. Our results strongly support that the accumulation of ECMs is essential for the development and progression of airway remodeling in asthma. Some components of ECMs may be novel therapeutic targets in the treatment of airway diseases by suppressing airway remodeling and inflammation.
It is also worth noting that the networks driven by transcription factors HIF-1α and FOXA1 and the signaling mediated by the interleukin-4 and 13 (IL4 and IL13) receptors were enriched for asthma (Figure 7). Previous studies revealed that HIF-1α was a master regulator of inflammation and was up-regulated in the immune cells and lung tissues of asthma patients [68,69,70].
In addition, hypoxia-induced HIF-1α could increase the expression of MUC5AC via NF-kB-mediated signaling pathway in asthma. Knocking down HIF-1α by siRNA decreased MUC5AC expression under hypoxia even in IL-1β-treated cells. The above results suggest that HIF-1α may be a therapeutic target for asthma [71].
Cytokine receptors binding to IL-4 activate the JAK-STAT signaling pathway that regulates the expression of downstream genes MMP1, ALOX15, CD36, VEGFA, FOS, NOS2, CLCA1, HIF-1α, TIMP1, and FN1 [72]. All target genes, except for HIF-1α, were identified in this study. Therefore, targeting this pathway through the inhibition of activating IL13 and IL4, and their receptors, and other pathway components should have therapeutic effects on asthma.

3. Materials and Methods

3.1. Microarray Gene Expression Data Acquisition

Publicly available microarray gene expression datasets for asthma were retrieved from the Gene Expression Omnibus Database (GEO) (http://www.ncbi.nlm.nih.gov/geo/) using the keyword “asthma”. The raw datasets were manually checked and only those that met the following criteria were included for subsequent analysis: 1) gene expression profiling in asthmatics and controls, 2) cell type: airway epithelial cell, but not nasal epithelium, 3) gene expression data were generated by a single-channel microarray platform (Affymetrix or Agilent chips), 4) availability of raw CEL or TXT files, and 5) samples with detailed descriptions.

3.2. Quality Control and Individual Microarray Dataset Preprocessing

Quality control (QC) analysis was performed for each microarray data using the arrayQualityMetrics package in R. Datasets were excluded from subsequent analysis when they failed to pass any one of the following assessments: 1) pairwise distances between arrays, 2) boxplots; 3) MA plots [73].
Briefly, the raw Affymetrix CEL files were normalized using the robust multi-array averaging (RMA) approach in the Affy package in Bioconductor Project with the following parameters: background correction, normalization and summarization [74], and returning log2 transformed intensities. Raw Agilent microarray data were normalized using the limma package, following background correction using the normexp method, quantile normalization and log2 transformation [75]. All analyses were performed with Bioconductor and R [76].

3.3. Data Integration and Batch Effect Removal

In this study, each dataset from all studies was preprocessed separately and then combined into one large dataset for further analysis following the steps illustrated in Figure 1. Probe identifiers from different microarray data were converted into Entrez gene IDs. If more than one probe mapped to a gene, the probe with the largest interquartile range (IQR) was selected as described by Letellier et al. [77], while a probe that mapped to multiple genes was excluded from further analysis.
The direct data integration (DDI) strategy presented in our previous publication [78] was applied to combine the multiple datasets in this study. Specifically, the shared genes among all microarray platforms used in this study were extracted, and then the multiple normalized datasets were merged into a large dataset (i.e., gene expression matrix) for shared genes.
To reduce potential study-specific batch effects, the preprocessed and normalized individual dataset was subjected to correction using the ComBat algorithm, which used the empirical Bayes method to adjust the extreme expression ratios and stabilized gene variances across all other genes and protected their reference from artifacts in the data [79]. The results of data integration and batch effect removal were visualized using clustering dendrograms, created by R and ggtree package [80].

3.4. Identification of Differentially Expressed Genes

Differentially expressed genes (DEGs) in integrated datasets were identified using limma package [75]. The linear fit, empirical Bayes (eBayes) statistics, and a false discovery rate (FDR) correction for all data were conducted using lmFit function. Genes with fold change > 1.2 and FDR < 0.05 were considered as DEGs between asthmatics and healthy controls. In addition, fold change > 1.5 and FDR < 0.05 were provided as well, hence we could identify genes which were more closely related to asthma.

3.5. Gene Ontology and Pathway Enrichment Analysis

To study the function of DEGs, enrichment of gene ontology (GO) was performed using DAVID (https://david.ncifcrf.gov/). The gene list was uploaded and analyzed using the annotation clustering for biological processes (BP), cellular components (CC), and molecular functions (MF). GO terms were considered to be significant when adjusted p < 0.01. The enrichment results were visualized using GOplot package [81]. Additionally, the enrichment analysis of Kyoto Encyclopedia of Genes and Genomes (KEGG)-based pathway was carried out using the GSEA online tool (http://www.broadinstitute.org/gsea/msigdb/index.jsp) [82]; KEGG pathways were considered to be significant when FDR < 0.05.
To identify up- or down-regulated terms based on DEGs, the Z-score proposed by Walter was adopted and calculated using the following formula [81]:
Z - s c o r e = N u p N d o w n C o u n t
where Nup, Ndown represent the number of up- or down-regulated genes between asthmatics and healthy controls, respectively. The count is the number of DEGs belonging to this term.
To investigate the biological pathway most relevant to asthma, functional pathway enrichment analyses of DEGs were performed using ToppGene Suite [83].

3.6. Protein-Protein Interaction Network Construction and Community Detection

To explore potential protein-protein interaction associated with asthma, DEGs with fold change > 1.5 were mapped to the PPI data using NetworkAnalyst [41] and its built-in IMEx Interactome, as well as the STRING online database [84]. The PPI network was visualized using Cytoscape V3.5.1 [85].
The asthma-relevant hub-genes were screened using the node degrees calculated in Cytoscape, and the detection of communities in the PPI network was performed with the InfoMap algorithm [56].

3.7. Target Gene Prediction of Key Transcription Factors and Regulatory MicroRNAs

Analysis of key transcription factors and their putative target genes was carried out using the GSEA online tool [82]. The statistical significance was determined using hypergeometric distribution and followed by Benjamini-Hochberg multiple testing. Significantly enriched TFs and microRNAs were considered when FDR < 0.05.

3.8. Chromosome Position Effect on Gene Expression

To study whether chromosome position affected gene expression patterns, all DEGs were mapped and enriched on chromosomes using the GSEA online tool. The mapping results (FDR < 0.01) were visualized using the karyoploteR package [86].

3.9. Connection of DEGs and Small Chemical Molecules

The Connectivity Map was a database of gene expression patterns from cultured human cells treated with bioactive small chemical molecules, in which 6100 groups of small molecule interference experiments and 7056 corresponding gene-expression profiles were stored [87]. CMap analysis was helpful in identifying bioactive molecules resulting in similar or adverse gene expression patterns. DEGs with fold change > 1.5 were subjected to connective mapping analysis using the PharmacoGx package [88], and the enrichment score (ES) was calculated. A molecule was considered beneficial for asthma when ES < −0.5, while it was harmful when ES > 0.5.

4. Conclusions

To investigate the characteristics of gene expression profiles of bronchial epithelium from asthmatic patients, data integration and systematic bioinformatics approaches were applied in this study. A total of 78 significantly up-regulated and 75 down-regulated genes were identified. Further analyses revealed that 10 differentially expressed genes are possibly instrumental in asthma, including up-regulated CEACAM5, CLCA1, POSTN, CPA3, SERPINB2, KRT6A, CD44, MUC5AC, and down-regulated LTF and MUC5B. Especially, the present study suggests that CD44, KRT6A, LTF, SERPINB2, MUC5B, and CEACAM5 could serve as hub genes in asthma-relevant PPI network and could further act as the biomarkers of asthma. Furthermore, the cross-talking of the IL-4/13 signaling pathway and networks intermediated by transcription factor HIF-1α and FOXA1 play a crucial role in the pathogenesis of asthma. Drug candidate prediction showed that compounds catechin, lomefloxacin, and boldine may be potential drugs for the treatment of asthma caused by inflammation and allergy. Our results elucidated the possible pathogenesis mechanisms of bronchial asthma and also provided several targets for further investigation.

Supplementary Materials

Supplementary materials can be found at https://www.mdpi.com/1422-0067/20/16/4037/s1.

Author Contributions

B.L. conceived the idea and supervised the work. X.N. and J.W. collected the datasets and completed the data integrating, X.N., J.W., and Y.H. carried out further microarray data analysis. J.T., Y.L., and M.L. wrote the original scripts and prepared the figures and tables, X.N. and B.X. wrote the original manuscript, and B.L. and Y.H. performed write-review and editing of the manuscript. All authors reviewed and approved the final version of the manuscript.

Funding

This work was funded by the Scientific and Technological Research Program of Chongqing Municipal Education Commission (No. KJQN201800523), Chongqing Research Program of Basic Research and Frontier Technology (No. cstc2015jcyjBX0142, cstc2016jcyjA0375), and Scientific Research Staring Foundation of Chongqing Normal University (No. 17XLB017).

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

GOGene Ontology
KEGGKyoto Encyclopedia of Genes and Genomes
DEGsDifferentially Expressed Genes
GSEAGene Set Enrichment Analysis
PPIProtein-protein Interaction
CMapConnectivity Map
ESEnrichment Score

References

  1. Willis-Owen, S.A.G.; Cookson, W.O.C.; Moffatt, M.F. The Genetics and Genomics of Asthma. Annu. Rev. Genomics Hum. Genet. 2018, 19, 223–246. [Google Scholar] [CrossRef] [PubMed]
  2. Lemmetyinen, R.; Karjalainen, J.; But, A.; Renkonen, R.; Pekkanen, J.; Toppila-Salmi, S.; Haukka, J. Higher mortality of adults with asthma: A 15-year follow-up of a population-based cohort. Allergy 2018, 73, 1479–1488. [Google Scholar] [CrossRef] [PubMed]
  3. Dhami, S.; Kakourou, A.; Asamoah, F.; Agache, I.; Lau, S.; Jutel, M.; Muraro, A.; Roberts, G.; Akdis, C.A.; Bonini, M. Allergen immunotherapy for allergic asthma: A systematic review and meta-analysis. Allergy 2017, 72, 1825–1848. [Google Scholar] [CrossRef] [PubMed]
  4. Huang, K.; Yang, T.; Xu, J.; Yang, L.; Zhao, J.; Zhang, X.; Bai, C.; Kang, J.; Ran, P.; Shen, H. Prevalence, risk factors, and management of asthma in China: A national cross-sectional study. Lancet 2019. [Google Scholar] [CrossRef]
  5. Forno, E.; Celedon, J.C. Epigenomics and Transcriptomics in the Prediction and Diagnosis of Childhood Asthma: Are We There Yet? Front. Pediatr. 2019, 7, 115. [Google Scholar] [CrossRef] [PubMed]
  6. Xiong, J.; Zhao, W.; Lin, Y.; Yao, L.; Huang, G.; Yu, C.; Dong, H.; Xiao, G.; Zhao, H.; Cai, S. Phosphorylation of low density lipoprotein receptor-related protein 6 is involved in receptor for advanced glycation end product-mediated β-catenin stabilization in a toluene diisocyanate-induced asthma model. Int. Immunopharmacol. 2018, 59, 187–196. [Google Scholar] [CrossRef]
  7. Lambrecht, B.N.; Hammad, H. The airway epithelium in asthma. Nat. Med. 2012, 18, 684. [Google Scholar] [CrossRef]
  8. Woodruff, P.G.; Boushey, H.A.; Dolganov, G.M.; Barker, C.S.; Yang, Y.H.; Donnelly, S.; Ellwanger, A.; Sidhu, S.S.; Dao-Pick, T.P.; Pantoja, C.; et al. Genome-wide profiling identifies epithelial cell genes associated with asthma and with treatment response to corticosteroids. Proc. Natl. Acad. Sci. USA 2007, 104, 15858–15863. [Google Scholar] [CrossRef] [Green Version]
  9. Kicic, A.; Hallstrand, T.S.; Sutanto, E.N.; Stevens, P.T.; Kobor, M.S.; Taplin, C.; Pare, P.D.; Beyer, R.P.; Stick, S.M.; Knight, D.A. Decreased fibronectin production significantly contributes to dysregulated repair of asthmatic epithelium. Am. J. Respir. Crit. Care Med. 2010, 181, 889–898. [Google Scholar] [CrossRef]
  10. Modena, B.D.; Tedrow, J.R.; Milosevic, J.; Bleecker, E.R.; Meyers, D.A.; Wu, W.; Bar-Joseph, Z.; Erzurum, S.C.; Gaston, B.M.; Busse, W.W.; et al. Gene expression in relation to exhaled nitric oxide identifies novel asthma phenotypes with unique biomolecular pathways. Am. J. Respir. Crit. Care Med. 2014, 190, 1363–1372. [Google Scholar] [CrossRef]
  11. Singhania, A.; Wallington, J.C.; Smith, C.G.; Horowitz, D.; Staples, K.J.; Howarth, P.H.; Gadola, S.D.; Djukanović, R.; Woelk, C.H.; Hinks, T.S. Multitissue transcriptomics delineates the diversity of airway T cell functions in asthma. Am. J. Respir. Cell Mol. Biol. 2018, 58, 261–270. [Google Scholar] [CrossRef] [PubMed]
  12. Karczewski, K.J.; Snyder, M.P. Integrative omics for health and disease. Nat. Rev. Genet. 2018, 19, 299. [Google Scholar] [CrossRef] [PubMed]
  13. Rung, J.; Brazma, A. Reuse of public genome-wide gene expression data. Nat. Rev. Genet. 2013, 14, 89. [Google Scholar] [CrossRef] [PubMed]
  14. Wagener, A.H.; Zwinderman, A.H.; Luiten, S.; Fokkens, W.J.; Bel, E.H.; Sterk, P.J.; van Drunen, C.M. The impact of allergic rhinitis and asthma on human nasal and bronchial epithelial gene expression. PLoS ONE 2013, 8, e80257. [Google Scholar] [CrossRef] [PubMed]
  15. Singhania, A.; Rupani, H.; Jayasekera, N.; Lumb, S.; Hales, P.; Gozzard, N.; Davies, D.E.; Woelk, C.H.; Howarth, P.H. Altered Epithelial Gene Expression in Peripheral Airways of Severe Asthma. PLoS ONE 2017, 12, e0168680. [Google Scholar] [CrossRef] [PubMed]
  16. Christenson, S.A.; Steiling, K.; van den Berge, M.; Hijazi, K.; Hiemstra, P.S.; Postma, D.S.; Lenburg, M.E.; Spira, A.; Woodruff, P.G. Asthma-COPD overlap. Clinical relevance of genomic signatures of type 2 inflammation in chronic obstructive pulmonary disease. Am. J. Respir. Crit. Care Med. 2015, 191, 758–766. [Google Scholar] [CrossRef]
  17. Yang, I.V.; Richards, A.; Davidson, E.J.; Stevens, A.D.; Kolakowski, C.A.; Martin, R.J.; Schwartz, D.A. The Nasal Methylome: A Key to Understanding Allergic Asthma. Am. J. Respir. Crit. Care Med. 2017, 195, 829–831. [Google Scholar] [CrossRef] [Green Version]
  18. Bindea, G.; Mlecnik, B.; Hackl, H.; Charoentong, P.; Tosolini, M.; Kirilovsky, A.; Fridman, W.-H.; Pagès, F.; Trajanoski, Z.; Galon, J. ClueGO: A Cytoscape plug-in to decipher functionally grouped gene ontology and pathway annotation networks. Bioinformatics 2009, 25, 1091–1093. [Google Scholar] [CrossRef]
  19. Greer, E.L.; Brunet, A. FOXO transcription factors at the interface between longevity and tumor suppression. Oncogene 2005, 24, 7410–7425. [Google Scholar] [CrossRef] [Green Version]
  20. Michnick, S.W. The connectivity map. Nat. Chem. Biol. 2006, 2, 663–664. [Google Scholar] [CrossRef]
  21. Brooks, W.H.; Renaudineau, Y. Epigenetics and autoimmune diseases: The X chromosome-nucleolus nexus. Front. Genet. 2015, 6, 22. [Google Scholar] [CrossRef]
  22. Rava, M.; Czachorowski, M.J.; Silverman, D.; Márquez, M.; Kishore, S.; Tardón, A.; Serra, C.; García-Closas, M.; Garcia-Closas, R.; Carrato, A. Asthma status is associated with decreased risk of aggressive urothelial bladder cancer. Int. J. Cancer. 2018, 142, 470–476. [Google Scholar] [CrossRef]
  23. Fish, J.E.; Ankin, M.G.; Adkinson, N.F., Jr.; Peterman, V.I. Indomethacin modification of immediate-type immunologic airway responses in allergic asthmatic and non-asthmatic subjects: Evidence for altered arachidonic acid metabolism in asthma. Am. Rev. Respir. Dis. 1981, 123, 609–614. [Google Scholar] [PubMed]
  24. Cottrell, L.; Neal, W.A.; Ice, C.; Perez, M.K.; Piedimonte, G. Metabolic abnormalities in children with asthma. Am. J. Respir. Crit. Care Med. 2011, 183, 441–448. [Google Scholar] [CrossRef] [PubMed]
  25. Woods, R.; Raven, J.; Walters, E.; Abramson, M.; Thien, F. Fatty acid levels and risk of asthma in young adults. Thorax 2004, 59, 105–110. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  26. Devereux, G.; Seaton, A. Diet as a risk factor for atopy and asthma. J. Allergy Clin. Immunol. 2005, 115, 1109–1117. [Google Scholar] [CrossRef] [PubMed]
  27. Mahn, K.; Ojo, O.O.; Chadwick, G.; Aaronson, P.I.; Ward, J.P.; Lee, T.H. Ca2+ homeostasis and structural and functional remodelling of airway smooth muscle in asthma. Thorax 2010, 65, 547–552. [Google Scholar] [CrossRef] [PubMed]
  28. Fritscher, L.; Chapman, K.R. Seretide: A pharmacoeconomic analysis. J. Med. Econ. 2008, 11, 555–570. [Google Scholar] [CrossRef] [PubMed]
  29. Paik, J.-h.; Ding, Z.; Narurkar, R.; Ramkissoon, S.; Muller, F.; Kamoun, W.S.; Chae, S.-S.; Zheng, H.; Ying, H.; Mahoney, J. FoxOs cooperatively regulate diverse pathways governing neural stem cell homeostasis. Cell stem cell 2009, 5, 540–553. [Google Scholar] [CrossRef]
  30. Schroeder, B.W.; Verhaeghe, C.; Park, S.-W.; Nguyenvu, L.T.; Huang, X.; Zhen, G.; Erle, D.J. AGR2 is induced in asthma and promotes allergen-induced mucin overproduction. Am. J. Respir. Cell Mol. Biol. 2012, 47, 178–185. [Google Scholar] [CrossRef]
  31. Shi, N.; Zhang, J.; Chen, S.-Y. Runx2, a novel regulator for goblet cell differentiation and asthma development. FASEB J. 2016, 31, 412–420. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  32. Makrinioti, H.; Toussaint, M.; Jackson, D.J.; Walton, R.P.; Johnston, S.L. Role of interleukin 33 in respiratory allergy and asthma. Lancet Respir. Med. 2014, 2, 226–237. [Google Scholar] [CrossRef]
  33. Zhou, X.; Kinlough, C.L.; Hughey, R.P.; Jin, M.; Inoue, H.; Etling, E.; Modena, B.D.; Kaminski, N.; Bleecker, E.R.; Meyers, D.A. Sialylation of MUC4β N-glycans by ST6GAL1 orchestrates human airway epithelial cell differentiation associated with type-2 inflammation. JCI insight 2019, 4. [Google Scholar] [CrossRef] [PubMed]
  34. Holgate, S.T.; Holloway, J.; Wilson, S.; Bucchieri, F.; Puddicombe, S.; Davies, D.E. Epithelial–mesenchymal communication in the pathogenesis of chronic asthma. Proc. Am. Thorac. Soc. 2004, 1, 93–98. [Google Scholar] [CrossRef]
  35. Feng, M.-J.; Shi, F.; Qiu, C.; Peng, W.-K. MicroRNA-181a,-146a and-146b in spleen CD4+ T lymphocytes play proinflammatory roles in a murine model of asthma. Int. Immunopharmacol. 2012, 13, 347–353. [Google Scholar] [CrossRef] [PubMed]
  36. Mohamed, J.S.; Lopez, M.A.; Boriek, A.M. Mechanical stretch up-regulates microRNA-26a and induces human airway smooth muscle hypertrophy by suppressing glycogen synthase kinase-3β. J. Biol. Chem. 2010, 285, 29336–29347. [Google Scholar] [CrossRef] [PubMed]
  37. Thomas, G.; Jacobs, K.B.; Yeager, M.; Kraft, P.; Wacholder, S.; Orr, N.; Yu, K.; Chatterjee, N.; Welch, R.; Hutchinson, A. Multiple loci identified in a genome-wide association study of prostate cancer. Nat. Genet. 2008, 40, 310. [Google Scholar] [CrossRef]
  38. Thomas, G.; Jacobs, K.B.; Kraft, P.; Yeager, M.; Wacholder, S.; Cox, D.G.; Hankinson, S.E.; Hutchinson, A.; Wang, Z.; Yu, K. A multistage genome-wide association study in breast cancer identifies two new risk alleles at 1p11. 2 and 14q24. 1 (RAD51L1). Nat. Genet. 2009, 41, 579. [Google Scholar] [CrossRef]
  39. Celedon, J.C.; Soto-Quiros, M.E.; Avila, L.; Lake, S.L.; Liang, C.; Fournier, E.; Spesny, M.; Hersh, C.P.; Sylvia, J.S.; Hudson, T.J.; et al. Significant linkage to airway responsiveness on chromosome 12q24 in families of children with asthma in Costa Rica. Hum. Genet. 2007, 120, 691–699. [Google Scholar] [CrossRef]
  40. Ferreira, M.A.; Matheson, M.C.; Duffy, D.L.; Marks, G.B.; Hui, J.; Le Souëf, P.; Danoy, P.; Baltic, S.; Nyholt, D.R.; Jenkins, M. Identification of IL6R and chromosome 11q13. 5 as risk loci for asthma. Lancet 2011, 378, 1006–1014. [Google Scholar] [CrossRef]
  41. Zhou, G.; Soufan, O.; Ewald, J.; Hancock, R.E.; Basu, N.; Xia, J. NetworkAnalyst 3.0: A visual analytics platform for comprehensive gene expression profiling and meta-analysis. Nucleic Acids Res. 2019. [Google Scholar] [CrossRef]
  42. Peroni, D.G.; Djukanovic, R.; Bradding, P.; Feather, I.H.; Montefort, S.; Howarth, P.H.; Jones, D.B.; Holgate, S.T. Expression of CD44 and integrins in bronchial mucosa of normal and mildly asthmatic subjects. Eur. Respir. J. 1996, 9, 2236–2242. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  43. Sano, K.; Yamauchi, K.; Hoshi, H.; Honma, M.; Tamura, G.; Shirato, K. CD44 expression on blood eosinophils is a novel marker of bronchial asthma. Int. Arch. Allergy Immunol. 1997, 114, 67–71. [Google Scholar] [CrossRef] [PubMed]
  44. Lv, Y.; Yang, S.; Zhang, Z.; Cui, Y.; Quan, C.; Zhou, F.; Fang, Q.; Du, W.; Zhang, F.; Chang, J. Novel and recurrent keratin 6A (KRT6A) mutations in Chinese patients with pachyonychia congenita type 1. Br. J. Dermatol. 2009, 160, 1327–1329. [Google Scholar] [CrossRef] [PubMed]
  45. Hu, J.; Zhang, L.-C.; Song, X.; Lu, J.-R.; Jin, Z. KRT6 interacting with notch1 contributes to progression of renal cell carcinoma, and aliskiren inhibits renal carcinoma cell lines proliferation in vitro. Int. J. Clin. Exp. Pathol. 2015, 8, 9182. [Google Scholar] [PubMed]
  46. Bu, W.; Chen, J.; Morrison, G.D.; Huang, S.; Creighton, C.J.; Huang, J.; Chamness, G.C.; Hilsenbeck, S.G.; Roop, D.R.; Leavitt, A.D. Keratin 6a marks mammary bipotential progenitor cells that can give rise to a unique tumor model resembling human normal-like breast cancer. Oncogene 2011, 30, 4399. [Google Scholar] [CrossRef]
  47. Gao, P.; Grigoryev, D.; Breslin, L.; Cheadle, C.; Mathias, R.; Beaty, T.; Togias, A.; Barnes, K. Keratins: Important candidate genes for asthma and immune responsiveness to cockroach. J. Allergy Clin. Immunol. 2007, 119, S176. [Google Scholar] [CrossRef]
  48. Zhang, C.; Abudula, A.; Awuti, M.; Wang, H.; Aihemaiti, X.; Tusung, T.; Sulaiman, X.; Upur, H. Plasma proteins as potential targets of abnormal Savda syndrome in asthma patients treated with unique Uighur prescription. Exp. Ther. Med. 2017, 14, 267–275. [Google Scholar] [CrossRef]
  49. Fang, F.; Pan, J.; Li, Y.; Li, Y.; Feng, X.; Wang, J. Identification of potential transcriptomic markers in developing asthma: An integrative analysis of gene expression profiles. Mol. Immunol. 2017, 92, 38–44. [Google Scholar] [CrossRef]
  50. Kruzel, M.L.; Bacsi, A.; Choudhury, B.; Sur, S.; Boldogh, I. Lactoferrin decreases pollen antigen-induced allergic airway inflammation in a murine model of asthma. Immunology 2006, 119, 159–166. [Google Scholar] [CrossRef]
  51. Tanner, A.; Ward, J.; Wilson, S.; Howarth, P. S30 CEACAM5 (CD66e) mucosal immunoreactivity and its relationship to asthma. Thorax 2018, 73, A18. [Google Scholar]
  52. Adcock, I.M.; Loza, M.; Djukanovic, R.; Rowe, A.; Rao, N.; Chung, K.F.; Baribaud, F. Differential gene expression profiles of endobronchial biopsies from severe and non-severe asthmatics obtained in UBIOPRED. Am. J. Respir. Crit. Care Med. 2014, A2422. [Google Scholar]
  53. de Jong, E.; Bosco, A. Dissecting Asthma Transcriptomics: Does Site Matter? Am. J. Respir. Cell Mol. Biol. 2018, 58, 144–146. [Google Scholar] [CrossRef] [PubMed]
  54. Naumov, D.; Gassan, D.; Kotova, O.; Prikhodko, A.; Pinegina, E.; Perelman, J.; Kolosov, V.; Zhou, X.; Li, Q. Cold air alters MUC5AC and MUC5B expression in the airways of asthma patients. Eur. Respir. J. 2018, 52 (Suppl. 62), PA1272. [Google Scholar]
  55. Shrine, N.; Portelli, M.A.; John, C.; Artigas, M.S.; Bennett, N.; Hall, R.; Lewis, J.; Henry, A.P.; Billington, C.K.; Ahmad, A. Moderate-to-severe asthma in individuals of European ancestry: A genome-wide association study. Lancet Respir. Med. 2019, 7, 20–34. [Google Scholar] [CrossRef]
  56. Zaki, N.; Mora, A. A comparative analysis of computational approaches and algorithms for protein subcomplex identification. Sci. Rep. 2014, 4, 4262. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  57. Nagata, H.; Takekoshi, S.; Takagi, T.; Honma, T.; Watanabe, K. Antioxidative action of flavonoids, quercetin and catechin, mediated by the activation of glutathione peroxidase. Tokai J. Exp. Clin. Med. 1999, 24, 1–11. [Google Scholar]
  58. Yamane, T.; Nakatani, H.; Kikuoka, N.; Matsumoto, H.; Iwata, Y.; Kitao, Y.; Oya, K.; Takahashi, T. Inhibitory effects and toxicity of green tea polyphenols for gastrointestinal carcinogenesis. Cancer 1996, 77, 1662–1667. [Google Scholar] [CrossRef]
  59. Fujimura, Y.; Tachibana, H.; Maeda-Yamamoto, M.; Miyase, T.; Sano, M.; Yamada, K. Antiallergic tea catechin,(−)-epigallocatechin-3-O-(3-O-methyl)-gallate, suppresses FcεRI expression in human basophilic KU812 cells. J. Agric. Food Chem. 2002, 50, 5729–5734. [Google Scholar] [CrossRef]
  60. Kim, S.H.; Park, H.J.; Lee, C.M.; Choi, I.W.; Moon, D.O.; Roh, H.J.; Lee, H.K.; Park, Y.M. Epigallocatechin-3-gallate protects toluene diisocyanate-induced airway inflammation in a murine model of asthma. FEBS Lett. 2006, 580, 1883–1890. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  61. Patel, S.; Patel, V. Inhibitory effects of catechin isolated from Acacia catechu on ovalbumin induced allergic asthma model: Role of histidine decarboxylase. Nutr. & Food Sci. 2019, 49, 18–31. [Google Scholar]
  62. Grassi, C.; Albera, C.; Pozzi, E. Lomefloxacin versus amoxicillin in the treatment of acute exacerbations of chronic bronchitis: An Italian multicenter study. Am. J. Med. 1992, 92, S103–S107. [Google Scholar] [CrossRef]
  63. Todd, I.; Negm, O.H.; Reps, J.; Radford, P.; Figueredo, G.; McDermott, E.M.; Drewe, E.; Powell, R.J.; Bainbridge, S.; Hamed, M. A signalome screening approach in the autoinflammatory disease TNF receptor associated periodic syndrome (TRAPS) highlights the anti-inflammatory properties of drugs for repurposing. Pharmacol. Res. 2017, 125, 188–200. [Google Scholar] [CrossRef]
  64. de Lima, N.M.R.; Ferreira, E.d.O.; Fernandes, M.Y.S.; Lima, F.A.V.; Neves, K.R.T.; do Carmo, M.R.S.; de Andrade, G.M. Neuroinflammatory response to experimental stroke is inhibited by boldine. Behav. pharmacol. 2017, 28, 223–237. [Google Scholar] [CrossRef]
  65. Fernández, J.; Lagos, P.; Rivera, P.; Zamorano-Ponce, E. Effect of boldo (Peumus boldus Molina) infusion on lipoperoxidation induced by cisplatin in mice liver. Phytother. Res. 2009, 23, 1024–1027. [Google Scholar] [CrossRef]
  66. Liu, G.; Cooley, M.A.; Nair, P.M.; Donovan, C.; Hsu, A.C.; Jarnicki, A.G.; Haw, T.J.; Hansbro, N.G.; Ge, Q.; Brown, A.C. Airway remodelling and inflammation in asthma are dependent on the extracellular matrix protein fibulin-1c. J. Pathol. 2017, 243, 510–523. [Google Scholar] [CrossRef]
  67. Agache, I.; Akdis, C.; Jutel, M.; Virchow, J. Untangling asthma phenotypes and endotypes. Allergy 2012, 67, 835–846. [Google Scholar] [CrossRef]
  68. Ma, L.; Javier, C.; Al Kolla, R.; Moshensky, A.; Bojanowski, C.; Shin, J.; Crotty Alexander, L.; Chun, V. HIFα Subunits Alter Asthma Phenotype in a Murine Model. Am. J. Respir. Crit. Care Med. 2018, 197, A1332. [Google Scholar]
  69. Dewitz, C.; McEachern, E.; Shin, S.; Akong, K.; Nagle, D.G.; Broide, D.H.; Akuthota, P.; Alexander, L.E.C. Hypoxia-inducible factor-1α inhibition modulates airway hyperresponsiveness and nitric oxide levels in a BALB/c mouse model of asthma. Clin. Immunol. 2017, 176, 94–99. [Google Scholar] [CrossRef]
  70. Ahmad, T.; Kumar, M.; Mabalirajan, U.; Pattnaik, B.; Aggarwal, S.; Singh, R.; Singh, S.; Mukerji, M.; Ghosh, B.; Agrawal, A. Hypoxia response in asthma: Differential modulation on inflammation and epithelial injury. Am. J. Respir. Cell Mol. Biol. 2012, 47, 1–10. [Google Scholar] [CrossRef]
  71. Wu, S.; Li, H.; Yu, L.; Wang, N.; Li, X.; Chen, W. IL-1β upregulates Muc5ac expression via NF-κB-induced HIF-1α in asthma. Immunol. lett. 2017, 192, 20–26. [Google Scholar] [CrossRef] [PubMed]
  72. O’Shea, J.J.; Plenge, R. JAK and STAT signaling molecules in immunoregulation and immune-mediated disease. Immunity 2012, 36, 542–550. [Google Scholar] [CrossRef] [PubMed]
  73. Kauffmann, A.; Gentleman, R.; Huber, W. arrayQualityMetrics--a bioconductor package for quality assessment of microarray data. Bioinformatics 2009, 25, 415–416. [Google Scholar] [CrossRef] [PubMed]
  74. Wu, Z.; Irizarry, R.A. Preprocessing of oligonucleotide array data. Nat. Biotechnol. 2004, 22, 656–658, 656–658; author reply 658. [Google Scholar] [CrossRef] [PubMed]
  75. Ritchie, M.E.; Phipson, B.; Wu, D.; Hu, Y.; Law, C.W.; Shi, W.; Smyth, G.K. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015, 43, e47. [Google Scholar] [CrossRef] [PubMed]
  76. Huber, W.; Carey, V.J.; Gentleman, R.; Anders, S.; Carlson, M.; Carvalho, B.S.; Bravo, H.C.; Davis, S.; Gatto, L.; Girke, T.; et al. Orchestrating high-throughput genomic analysis with Bioconductor. Nat. Methods 2015, 12, 115–121. [Google Scholar] [CrossRef]
  77. Letellier, E.; Schmitz, M.; Ginolhac, A.; Rodriguez, F.; Ullmann, P.; Qureshi-Baig, K.; Frasquilho, S.; Antunes, L.; Haan, S. Loss of Myosin Vb in colorectal cancer is a strong prognostic factor for disease recurrence. Br. J. Cancer 2017, 117, 1689. [Google Scholar] [CrossRef]
  78. Yang, Q.; Wang, Y.; Zhang, S.; Tang, J.; Li, F.; Yin, J.; Li, Y.; Fu, J.; Li, B.; Luo, Y. Biomarker discovery for immunotherapy of pituitary adenomas: Enhanced robustness and prediction ability by modern computational tools. Int. J. Mol. Sci. 2019, 20, 151. [Google Scholar] [CrossRef]
  79. Leek, J.T.; Johnson, W.E.; Parker, H.S.; Jaffe, A.E.; Storey, J.D. The sva package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinformatics 2012, 28, 882–883. [Google Scholar] [CrossRef]
  80. Yu, G.; Smith, D.K.; Zhu, H.; Guan, Y.; Lam, T.T.Y. ggtree: An R package for visualization and annotation of phylogenetic trees with their covariates and other associated data. Methods Ecol. Evol. 2017, 8, 28–36. [Google Scholar] [CrossRef]
  81. Walter, W.; Sánchez-Cabo, F.; Ricote, M. GOplot: An R package for visually combining expression data with functional analysis. Bioinformatics 2015, 31, 2912–2914. [Google Scholar] [CrossRef] [PubMed]
  82. Reimand, J.; Isserlin, R.; Voisin, V.; Kucera, M.; Tannus-Lopes, C.; Rostamianfar, A.; Wadi, L.; Meyer, M.; Wong, J.; Xu, C. Pathway enrichment analysis and visualization of omics data using g: Profiler, GSEA, Cytoscape and EnrichmentMap. Nat. Protoc. 2019, 14, 482. [Google Scholar] [CrossRef] [PubMed]
  83. Chen, J.; Bardes, E.E.; Aronow, B.J.; Jegga, A.G. ToppGene Suite for gene list enrichment analysis and candidate gene prioritization. Nucleic Acids Res. 2009, 37, W305–311. [Google Scholar] [CrossRef] [PubMed]
  84. Szklarczyk, D.; Gable, A.L.; Lyon, D.; Junge, A.; Wyder, S.; Huerta-Cepas, J.; Simonovic, M.; Doncheva, N.T.; Morris, J.H.; Bork, P. STRING v11: Protein–protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res. 2018, 47, D607–D613. [Google Scholar] [CrossRef] [PubMed]
  85. Cline, M.S.; Smoot, M.; Cerami, E.; Kuchinsky, A.; Landys, N.; Workman, C.; Christmas, R.; Avila-Campilo, I.; Creech, M.; Gross, B. Integration of biological networks and gene expression data using Cytoscape. Nat. Protoc. 2007, 2, 2366. [Google Scholar] [CrossRef] [PubMed]
  86. Gel, B.; Serra, E. karyoploteR: An R/Bioconductor package to plot customizable genomes displaying arbitrary data. Bioinformatics 2017, 33, 3088–3090. [Google Scholar] [CrossRef] [PubMed]
  87. Lamb, J.; Crawford, E.D.; Peck, D.; Modell, J.W.; Blat, I.C.; Wrobel, M.J.; Lerner, J.; Brunet, J.P.; Subramanian, A.; Ross, K.N.; et al. The Connectivity Map: Using gene-expression signatures to connect small molecules, genes, and disease. Science 2006, 313, 1929–1935. [Google Scholar] [CrossRef] [PubMed]
  88. Smirnov, P.; Safikhani, Z.; El-Hachem, N.; Wang, D.; She, A.; Olsen, C.; Freeman, M.; Selby, H.; Gendoo, D.M.; Grossmann, P.; et al. PharmacoGx: An R package for analysis of large pharmacogenomic datasets. Bioinformatics 2016, 32, 1244–1246. [Google Scholar] [CrossRef] [PubMed]
Figure 1. The workflow of microarray data integration and subsequent analyses in this study. Quality control of each dataset was manually checked, and then 10 preprocessed and independent qualifying datasets were merged into a large dataset (designated as Merged DS) for further analysis. Batch effects were removed by using the algorithm of batch effect removal. Differentially expressed genes (DEGs) between the asthmatics and healthy controls were identified based on expression fold change > 1.2 and false discovery rate (FDR) < 0.05. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG)-based pathways were subsequently enriched and followed by protein-protein interaction (PPI) network construction, prediction of potential transcription factors (TFs) and microRNAs, chromosomal localization, and potential candidate chemical prediction for asthma treatment. The partial contents in this workflow were cited from published papers [18,19,20,21].
Figure 1. The workflow of microarray data integration and subsequent analyses in this study. Quality control of each dataset was manually checked, and then 10 preprocessed and independent qualifying datasets were merged into a large dataset (designated as Merged DS) for further analysis. Batch effects were removed by using the algorithm of batch effect removal. Differentially expressed genes (DEGs) between the asthmatics and healthy controls were identified based on expression fold change > 1.2 and false discovery rate (FDR) < 0.05. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG)-based pathways were subsequently enriched and followed by protein-protein interaction (PPI) network construction, prediction of potential transcription factors (TFs) and microRNAs, chromosomal localization, and potential candidate chemical prediction for asthma treatment. The partial contents in this workflow were cited from published papers [18,19,20,21].
Ijms 20 04037 g001
Figure 2. Exploration and visualization of the Merged dataset (abbr. DS) before and after batch correcting. Hierarchical clustering for samples was performed based on Euclidean distance and Complete linkage. Then, the results of clustering were visualized using clustering dendrograms in ggtree package. (A) Before batch effect removal, and (B) after batch effect removal.
Figure 2. Exploration and visualization of the Merged dataset (abbr. DS) before and after batch correcting. Hierarchical clustering for samples was performed based on Euclidean distance and Complete linkage. Then, the results of clustering were visualized using clustering dendrograms in ggtree package. (A) Before batch effect removal, and (B) after batch effect removal.
Ijms 20 04037 g002
Figure 3. Volcano plot for 153 differentially expressed genes between asthmatics and healthy controls. The up- and down-regulated genes are marked in red and green, respectively. Genes with fold change more than 1.5-fold are labeled by gene symbols.
Figure 3. Volcano plot for 153 differentially expressed genes between asthmatics and healthy controls. The up- and down-regulated genes are marked in red and green, respectively. Genes with fold change more than 1.5-fold are labeled by gene symbols.
Ijms 20 04037 g003
Figure 4. Bubble plot for GO term enrichment. The bubble size represents the number of DEGs in each GO term. Three levels of GO terms are the biological process (BP), cellular component (CC), and molecular function (MF).
Figure 4. Bubble plot for GO term enrichment. The bubble size represents the number of DEGs in each GO term. Three levels of GO terms are the biological process (BP), cellular component (CC), and molecular function (MF).
Ijms 20 04037 g004
Figure 5. Chromosome localization of nine DEGs potentially associated with asthma
Figure 5. Chromosome localization of nine DEGs potentially associated with asthma
Ijms 20 04037 g005
Figure 6. Protein-protein interaction network associated with asthma was constructed based on the DEGs with fold change > 1.5. The up- and down-regulated DEGs are marked in red and green circles, respectively. The solid lines and dashed lines represent the interactions predicted based on the IMEx interactome database and STRING database (v11.0), respectively. According to the degree of each node, six proteins with degree ≥ 10 were regarded as the hub proteins, and their corresponding genes were thought to be the hub genes, including CD44 (degree = 59), KRT6A (degree = 23), LTF (degree = 18), MUC5B (degree = 13), CEACAM5 (degree = 11), and SERPINB2 (degree = 10).
Figure 6. Protein-protein interaction network associated with asthma was constructed based on the DEGs with fold change > 1.5. The up- and down-regulated DEGs are marked in red and green circles, respectively. The solid lines and dashed lines represent the interactions predicted based on the IMEx interactome database and STRING database (v11.0), respectively. According to the degree of each node, six proteins with degree ≥ 10 were regarded as the hub proteins, and their corresponding genes were thought to be the hub genes, including CD44 (degree = 59), KRT6A (degree = 23), LTF (degree = 18), MUC5B (degree = 13), CEACAM5 (degree = 11), and SERPINB2 (degree = 10).
Ijms 20 04037 g006
Figure 7. Cross-talking pathways associated with asthma inferred by DEGs as well as literature.
Figure 7. Cross-talking pathways associated with asthma inferred by DEGs as well as literature.
Ijms 20 04037 g007
Table 1. Information summary of microarray datasets used in this study.
Table 1. Information summary of microarray datasets used in this study.
GEO IDMicroarray PlatformSample Size (A/H) 1Cell TypeAuthor (Year)
GSE470HG-U95Av212 (6/6)Epithelial cellsSpannhake W, et al. (2003)
GSE4302HG-U133_Plus_270 (42/28)Airway epithelial cellsWoodruff PG, et al. (2007) [8]
GSE18965HG-U133A16 (9/7)Bronchial epithelial cellsBeyer RP, et al. (2010) [9]
GSE41861HG-U133_Plus_281 (51/30)Bronchial epithelial cellsCheng DT, et al. (2015)
GSE44037HT_HG-U133_Plus_PM12 (6/6)Bronchial epitheliaWagener AH, et al. (2013) [14]
GSE63142GPL6480 (Agilent)155 (128/27)Bronchial epitheliaWenzel S, et al. (2014) [10]
GSE64913HG-U133_Plus_259 (22/37)Peripheral airway epitheliaSinghania A, et al. (2017) [15]
GSE67472HG-U133_Plus_2105 (62/43)Airway epitheliaChristenson SA, et al. (2015) [16]
GSE89809 HT_HG-U133_Plus_PM56 (38/18)Epithelial cellsSinghania A, et al. (2017) [11]
GSE104468 GPL21185 (Agilent)24 (12/12)Bronchial epitheliaRichards A, et al. (2017) [17]
1 A and H represent asthmatics and healthy controls, respectively.
Table 2. The significant DEGs with fold change > 1.5 between asthmatics and healthy controls
Table 2. The significant DEGs with fold change > 1.5 between asthmatics and healthy controls
GeneEntrez IDLog2(Fold Change)Asthmatics vs. Healthy ControlsFDR 1
CEACAM510481.13Up7.67 × 10−23
CLCA111791.58Up5.43 × 10−22
POSTN106311.33Up7.83 × 10−22
CPA313591.26Up1.28 × 10−21
SERPINB250551.14Up9.55 × 10−20
LTF4057−0.75Down4.60 × 10−17
MUC5B727897−0.89Down2.10 × 10−13
KRT6A38530.61Up4.89 × 10−12
CD449600.60Up9.03 × 10−9
MUC5AC45860.59Up1.03 × 10−6
1 FDR (false discovery rate) refers to the BH-adjusted p-value returned by eBayes function in limma package.
Table 3. GO (gene ontology) term enrichment of DEGs between the asthmatics and healthy controls.
Table 3. GO (gene ontology) term enrichment of DEGs between the asthmatics and healthy controls.
GO TermDescriptionCount 1Z-scoreAdj. p 2Category
GO:0042221Response to chemical stimulus34−0.342.64 × 10−4BP
GO:0032501Multicellular organismal process690.121.74 × 10−3BP
GO:0018149Peptide cross-linking60.002.50 × 10−3BP
GO:0005576Extracellular region510.421.67 × 10−8CC
GO:0044421Extracellular region part310.184.21 × 10−7CC
GO:0005615Extracellular space21−0.654.33 × 10−4CC
GO:0030141Secretory granule100.632.75 × 10−3CC
GO:0031012Extracellular matrix130.284.33 × 10−3CC
GO:0000267Cell fraction24−0.826.39 × 10−3CC
GO:0005624Membrane fraction20−0.896.10 × 10−3CC
GO:0005578Proteinaceous extracellular matrix120.005.90 × 10−3CC
GO:0005626Insoluble fraction20−0.897.52 × 10−3CC
1 Count denotes the number of DEGs in this GO term; 2 Adj. p refers to the adjusted p-values computed by DAVID.
Table 4. Enriched KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways for DEGs using Molecular Signatures Database from GSEA (Gene Set Enrichment Analysis) online tool.
Table 4. Enriched KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways for DEGs using Molecular Signatures Database from GSEA (Gene Set Enrichment Analysis) online tool.
KEGG PathwayFDR 1Expression Pattern 2Z-scoreGene
(1) Pathways in cancer3.24 × 10−53↑+ 7↓−1.26NOS2, MMP1, VEGFA, DAPK1, FOS, KIT, FN1, RUNX1T1, EPAS1, AR
(2) Arachidonic acid metabolism1.42 × 10−43↑+ 2↓0.45CYP2J2, ALOX15, GPX3, HPGDS, PTGS1
(3) Linoleic acid metabolism7.69 × 10−32↑+ 1↓0.58CYP2J2, ALOX15, AKR1B10
(4) Calcium signaling pathway1.17 × 10−21↑+ 4↓−1.34NOS2, ITPR1, AVPR1A, PTGFR, TRPC1
(5) Aldosterone-regulated sodium reabsorption1.17 × 10−20↑+ 3↓−1.73IRS2, INSR, SCNN1G
(6) Bladder cancer1.17 × 10−21↑+ 2↓−0.58MMP1, VEGFA, DAPK1
(7) Arginine and proline metabolism1.95 × 10−22↑+ 1↓0.58NOS2, ODC1, PYCR1
(8) PPAR signaling pathway3.34 × 10−22↑+ 1↓0.58MMP1, CD36, FABP6
(9) Leishmania infection3.39 × 10−21↑+ 2↓−0.58NOS2, FOS, HLA-DQB1
(10) Cytokine-cytokine receptor interaction3.51 × 10−22↑+ 3↓−0.45VEGFA, KIT, CXCL2, CXCL6, CSF2RB
(11) ECM-receptor interaction4.39 × 10−22↑+ 1↓0.58FN1, CD36, CD44
(12) Hematopoietic cell lineage4.62 × 10−23↑+ 0↓1.73KIT, CD36, CD44
1 FDR denotes the FDR q-value provided by GSEA online tool; 2 ↑ refers to the status of up-regulated DEGs, and ↓ refers to the status of down-regulated DEGs.
Table 5. List of small molecules enriched from CMap database.
Table 5. List of small molecules enriched from CMap database.
MoleculesEnrichment Scorep-Value
Catechin−0.79950.0066
Lomefloxacin−0.71010.0072
Boldine−0.66350.0250
Prestwick-10820.69440.0099
Ricinine0.72970.0102
Milrinone0.77510.0153
Econazole0.79810.0170
Acetohexamide0.65420.0181
Cefsulodin0.79000.0192
Nifuroxazide0.75670.0199
Alimemazine0.77440.0289
Progesterone0.72580.0309
Zoxazolamine0.72060.0390
Colistin0.65030.0404
Methapyrilene0.76180.0427
Tiapride0.68320.0436
Fluocinonide0.74510.0466
Ganciclovir0.73890.0466
Quinisocaine0.67900.0479
Mexiletine0.76150.0492
Table 6. The enriched pathways from BIOCARTA, REACTOME and Pathway Interaction Database.
Table 6. The enriched pathways from BIOCARTA, REACTOME and Pathway Interaction Database.
Pathway IDNameDatabaseFDRCount 1
M5889Genes encoding ECM and ECM-associated proteinsBIOCARTA (v6.0)1.40 × 10−226
M5885ECM-affiliated proteins, regulators and secreted factors BIOCARTA (v6.0)2.42 × 10−220
1470923Interleukin−4 and 13 signalingREACTOME1.40 × 10−28
138045HIF−1-alpha transcription factor networkPathway Interaction Database2.15 × 10−26
137979FOXA1 transcription factor networkPathway Interaction Database2.23 × 10−25
1268737Termination of O-glycan biosynthesisREACTOME2.61 × 10−24
1 Count denotes the number of DEGs in this term.

Share and Cite

MDPI and ACS Style

Nie, X.; Wei, J.; Hao, Y.; Tao, J.; Li, Y.; Liu, M.; Xu, B.; Li, B. Consistent Biomarkers and Related Pathogenesis Underlying Asthma Revealed by Systems Biology Approach. Int. J. Mol. Sci. 2019, 20, 4037. https://doi.org/10.3390/ijms20164037

AMA Style

Nie X, Wei J, Hao Y, Tao J, Li Y, Liu M, Xu B, Li B. Consistent Biomarkers and Related Pathogenesis Underlying Asthma Revealed by Systems Biology Approach. International Journal of Molecular Sciences. 2019; 20(16):4037. https://doi.org/10.3390/ijms20164037

Chicago/Turabian Style

Nie, Xiner, Jinyi Wei, Youjin Hao, Jingxin Tao, Yinghong Li, Mingwei Liu, Boying Xu, and Bo Li. 2019. "Consistent Biomarkers and Related Pathogenesis Underlying Asthma Revealed by Systems Biology Approach" International Journal of Molecular Sciences 20, no. 16: 4037. https://doi.org/10.3390/ijms20164037

APA Style

Nie, X., Wei, J., Hao, Y., Tao, J., Li, Y., Liu, M., Xu, B., & Li, B. (2019). Consistent Biomarkers and Related Pathogenesis Underlying Asthma Revealed by Systems Biology Approach. International Journal of Molecular Sciences, 20(16), 4037. https://doi.org/10.3390/ijms20164037

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop