Microarray data analysis to identify differentially expressed genes and biological pathways associated with asthma.

The present study aimed to identify differentially expressed genes (DEGs) and biological processes (BPs) associated with asthma. DEGs between allergic asthma and healthy controls were screened from GSE15823. Gene Ontology and Kyoto Encyclopedia of Genes and Genomes pathway enrichment analyses were performed using Database for Annotation, Visualization and Integrated Discovery. A protein-protein interaction network was constructed, followed by module mining and functional analysis. Additionally, GSE41649 was downloaded to validate the reliability of the results. In GSE41649, DEGs were identified and compared with key DEGs identified in GSE15823. A total of 43 upregulated and 275 downregulated DEGs were obtained from GSE15823. Upregulated DEGs, such as nitric oxide synthase 2 (NOS2), were enriched in BPs related to oxidation reduction. Downregulated DEGs, such as chemokine (C-C motif) ligand 19 (CCL21) and Cys-X-Cys ligand (CXCL9), were enriched in immune response-related BPs. Protein tyrosine phosphatase receptor type C (PTPRC), CCL21, and CXCL9 were identified as hub genes. The DEGs in module 1 were significantly involved in the chemokine signaling pathway (P<0.05). The expression of the key genes obtained in GSE15823 demonstrated the same variation directions in the two datasets. The immune response, oxidants and nitric oxide metabolic pathways may have important roles in the progression of asthma. DEGs of PTPRC, CCL21, CXCL9 and NOS2 may be the potential targets for asthma diagnosis and treatment.

Introduction

Asthma is a severe inflammatory disorder of the airways that is characterized by hyperresponsiveness, obstruction and inflammation of the airways (1). It is a serious global health problem and individuals of all ages may be affected by this chronic airway disorder (2). The prevalence of asthma is rising globally (from 183 million in 1990 to ~242 million in 2013), particularly among children in industrialized countries (3,4). However, there is no cure for asthma at present, although symptoms may be improved with certain medications (5). Therefore, an improved understanding of the molecular mechanisms of this disease is required to discover novel therapeutic targets for patients with asthma.

Asthma is believed to be caused by the interaction of environmental and genetic factors (1). It has been suggested that environmental exposure, such as smoking, and various home allergens, including mold, dust mites and pollen, are associated with asthma (6). Estimates of heritability indicate that 35–80% of the variation in risk of developing asthma is due to genetic variation (7). For example, glutathione S-transferase mu 1, interleukin 10 (IL-10), cytotoxic T-lymphocyte-associated protein 4, leukotriene C4 synthase, IL-4 receptor and ADAM metallopeptidase domain 33 have been implicated in asthma (8). T-box 21 has been suggested to have a vital role in asthma pathobiology via regulation of T-helper 1 lineage commitment and interferon-γ generation (9). Furthermore, intracellular signaling component-associated pathways have been implicated in asthma pathogenesis (10). Although some potential biomarkers and functional pathways associated with asthma have been studied, current knowledge is not adequate enough for the effective diagnosis and treatment of this disease.

Bioinformatics has a vital role in identifying the underlying genetic basis of human disease (11). A study by Laprise et al (10) compared the expression of genes and transcripts between healthy controls and allergic asthma subjects in 2004 and deposited the gene expression profile, GSE15823, in the Gene Expression Omnibus (ncbi.nlm.nih.gov/geo/). A study by Vaillancourt et al (12) verified several differentially expressed genes (DEGs) identified by Laprise et al (10) as well as by Chamberland et al (13). However, Gene Ontology (GO), pathway functional enrichment and protein-protein interaction (PPI) networks based on the DEGs were not analyzed.

In the present study, bioinformatics methods were employed to identify DEGs in lung tissue samples from patients with asthma based on the gene expression profile of GSE15823. In addition, GO functions, pathways and PPI networks were analyzed to explore the potential functions, pathways and key genes involved in asthma. Furthermore, module mining and functional analysis from PPI networks were conducted. Notably, data validation based on another gene expression profile dataset was performed in order to validate the findings of the present study.

Materials and methods

Analysis of Affymetrix microarray data

The gene expression profile of GSE15823, based on the GPL8300 platform of Affymetrix Human Genome U95 Version 2 Array, was downloaded from the National Center of Biotechnology Information GEO database (ncbi.nlm.nih.gov/geo/). According to the criteria of the American Thoracic Society for the diagnosis of asthma (14), A total of 4 healthy subjects without a history of allergy or asthma and 4 patients with asthma meeting the criteria were recruited to the study by Laprise et al (10). Lung tissue samples were obtained from bronchial biopsies of the 8 subjects. In the present study, gene expression data of the 8 lung tissue samples were used for further analysis.

Data preprocessing and identification of DEGs

Raw Affymetrix data (in the form of CEL files) were downloaded and converted into expression values, followed by background correction. Subsequently, gene expression data were normalized using the robust multiarray average algorithm (15) of the Affy package (16) in Bioconductor software (version 3.4; http://www.bioconductor.org/packages/release/bioc/), followed by probe summarization.

Student's t-tests (17) were used to identify DEGs between asthma and healthy control samples. Adjusted P-values were obtained after the t-tests were performed and the fold-change (FC) was also calculated. An adjusted P-value of P<0.1 and |log2FC|>0.5 were selected as the criteria for DEG screening.

GO and pathway functional enrichment analyses of DEGs

Database for Annotation, Visualization and Integrated Discovery (DAVID; david.ncifcrf.gov) (18) is a tool that provides an extensive set of functional annotation. GO database (geneontology.org) (19) is a tool used for functional unification of large-scale genomics and includes three categories: Biological process (BP), molecular function and cellular component. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database (20) covers various biochemical pathways. In the present study, GO BP functional and KEGG pathway enrichment analysis for upregulated and downregulated DEGs were performed with DAVID. The modified Fisher's exact P-value of P<0.1 and count >2 were used as thresholds.

PPI network construction and module functional analysis

The Search Tool for the Retrieval of Interacting Genes (STRING) (21) database provides experimental and predicted interaction information. In the present study, the STRING online tool was used to analyze the PPI pairs of DEGs with a combined score of >0.4. The PPI network was then visualized by Cytoscape (version 3.2.0) (22) and the proteins (genes) with the highest connectivity degrees in the PPI network were regarded as hub proteins (genes).

Modules may be defined in a data-driven fashion by identifying subgroups of genes sharing similar expression patterns across multiple conditions (23). From the PPI network, the significant modules with the top three Molecular Complex Detection (MCODE) scores were extracted for analysis using MCODE of Cytoscape (24) with default thresholds. These thresholds included: Degree cutoff, 2; node score cutoff, 0.2; k-core, 2; and maximum depth, 100. Subsequently, the DEGs in modules underwent GO and pathway enrichment analyses.

Data validation

In order to validate the reliability of the results of the present study, expression profile data of GSE41649 (13) were downloaded from the GEO database, based on the platform of GPL96 (HG-U133A) Affymetrix Human Genome U133A (Affymetrix Inc., Santa Clara, CA, USA). Additionally, 8 samples, including 4 healthy controls without a history of allergy or asthma and 4 patients with allergic asthma, were included in the dataset. DEGs between allergic asthma and control samples were identified and compared with the key DEGs identified in GSE15823. The pathways and GO functions enriched by these DEGs were analyzed using DAVID. P<0.05 and a count of >2 were considered to indicate a statistically significant result.

Results

Identification of DEGs

As demonstrated in Fig. 1, the raw expression data were well-normalized after preprocessing. A total of 318 DEGs were obtained between asthma samples and healthy subjects after microarray analysis, including 43 upregulated and 275 downregulated DEGs.

Figure 1.

Box plots of the normalization of expression data. Data (A) prior to and (B) following normalization. Log intensity, gene expression values.

GO and KEGG pathway enrichment analysis

The enriched functions of upregulated and downregulated DEGs are demonstrated in Table I. The upregulated DEGs, such as nitric oxide synthase 2 inducible (NOS2), were enriched in the BPs associated with oxidation reduction and nitric oxide metabolism. The downregulated DEGs, such as chemokine (C-C motif) ligand 21 (CCL21) and Cys-X-Cys ligand (CXCL9), were predominantly enriched in the BPs associated with the immune response.

Table I.

GO functional enrichment of upregulated and downregulated DEGs.

Category	Description	Count	P-value
Upregulated DEGs
BP	GO: 0055114: Oxidation reduction	6	1.62E-02
BP	GO: 0043523: Regulation of neuron apoptosis	3	1.91E-02
BP	GO: 0006809: Nitric oxide biosynthetic process	2	2.80E-02
BP	GO: 0046209: Nitric oxide metabolic process	2	3.03E-02
BP	GO: 0042981: Regulation of apoptosis	6	3.88E-02
BP	GO: 0043067: Regulation of programmed cell death	6	4.03E-02
BP	GO: 0010941: Regulation of cell death	6	4.08E-02
BP	GO: 0031214 Biomineral formation	2	8.40E-02
Downregulated DEGs
BP	GO: 0006955: Immune response	44	1.89E-14
BP	GO: 0002449: Lymphocyte mediated immunity	11	1.86E-07
BP	GO: 0016064: Immunoglobulin mediated immune response	10	1.98E-07
BP	GO: 0019724: B cell mediated immunity	10	2.74E-07
BP	GO: 0002252: Immune effector process	14	2.87E-07
BP	GO: 0002460: Adaptive immune response based on somatic recombination of immune receptors built from immunoglobulin superfamily domains	11	4.66E-07
BP	GO: 0002250: Adaptive immune response	11	4.66E-07
BP	GO: 0006952: Defense response	29	8.73E-07
BP	GO: 0002443: Leukocyte mediated immunity	11	1.32E-06
BP	GO: 0002684: Positive regulation of immune system process	17	1.85E-06

GO, gene ontology; DEGs, differentially expressed genes; BP, biological process; Description, the enriched GO function of DEGs; Count, number of DEGs.

The KEGG pathways of upregulated and downregulated DEGs are demonstrated in Table II. Upregulated DEGs were enriched in caffeine metabolism. Downregulated DEGs were enriched in immune-associated pathways, such as systemic lupus erythematosus, allograft rejection, and complement and coagulation cascades.

Table II.

KEGG pathway enrichment of upregulated and downregulated DEGs.

Category	Description	Count	P-value
Upregulated DEGs
KEGG	hsa00232: Caffeine metabolism	2	1.91E-02
Downregulated DEGs
KEGG	hsa05322: Systemic lupus erythematosus	12	7.52E-06
KEGG	hsa04672: Intestinal immune network for IgA production	8	7.43E-05
KEGG	hsa05330: Allograft rejection	7	1.03E-04
KEGG	hsa05416: Viral myocarditis	9	1.27E-04
KEGG	hsa05332: Graft-versus-host disease	7	1.63E-04
KEGG	hsa04940: Type I diabetes mellitus	7	2.49E-04
KEGG	hsa05310: Asthma	6	3.39E-04
KEGG	hsa04610: Complement and coagulation cascades	8	6.54E-04
KEGG	hsa05320: Autoimmune thyroid disease	7	7.29E-04
KEGG	hsa04612: Antigen processing and presentation	8	1.96E-03
KEGG	hsa04514: Cell adhesion molecules	10	2.06E-03
KEGG	hsa05340: Primary immunodeficiency	5	6.48E-03
KEGG	hsa00071: Fatty acid metabolism	5	1.04E-02
KEGG	hsa04060: Cytokine-cytokine receptor interaction	13	1.05E-02
KEGG	hsa04062: Chemokine signaling pathway	10	1.92E-02
KEGG	hsa04640: Hematopoietic cell lineage	6	3.77E-02
KEGG	hsa05020: Prion diseases	4	3.91E-02

KEGG, Kyoto Encyclopedia of Genes and Genomes; DEGs, differentially expressed genes; Description, the enriched GO function of DEGs; Count, number of DEGs.

PPI network construction

In total, 220 DEGs with 381 PPI interaction pairs were identified in the PPI network (Fig. 2). The top 20 nodes with the highest degrees in the PPI networks are demonstrated in Table III and were defined as hub proteins (genes). These hub proteins (genes) included cluster of differentiation 8a molecule (CD8A; degree, 30), protein tyrosine phosphatase receptor type C (PTPRC; degree, 26), fibronectin 1 (FN1; degree, 14), CCL21 (degree, 13) and CXCL9 (degree, 10).

Figure 2.

Protein-protein interaction network of DEGs. Red node represents upregulated DEGs; green node represents downregulated DEGs; the different shade of color represents the significance of the difference, with deeper color indicating a more significant difference. DEGs, differentially expressed genes.

Table III.

Top 20 genes with the highest connectivity degrees in the protein-protein interaction network.

Gene name	Degree	Adjusted P-value	Log2FC
CD8A	30	0.23	−0.66
PTPRC	26	0.22	−0.50
FN1	14	0.24	−0.56
CCL19	13	0.26	−0.76
CCL5	12	0.24	−0.78
CCL21	11	0.20	−0.60
CXCR4	11	0.20	−1.04
CXCL9	10	0.32	−0.86
ADCY7	10	0.17	−0.80
CXCR6	10	0.25	−0.52
CD2	10	0.20	−0.60
CXCL6	9	0.20	−1.07
ITK	8	0.23	−0.53
CD247	8	0.23	−0.55
PENK	8	0.22	−0.65
H2AFV	8	0.21	−0.56
POSTN	7	0.21	1.81
SERPING1	7	0.21	−0.72
KLRB1	7	0.23	−0.74
COL3A1	7	0.29	−0.61

Degree, connectivity degree; FC, fold-change.

Module mining and functional analysis

Module analysis demonstrated that there were 9, 6 and 6 nodes in the 3 modules, respectively, and the DEGs involved in these modules were all downregulated (Fig. 3). Their MCODE scores were 9, 6 and 8, respectively.

Figure 3.

Modules in the protein-protein interaction network. There are (A) 9 nodes in module 1, (B) 6 nodes in module 2 and (C) 6 nodes in module 3.

As demonstrated in Table IV, DEGs in module 1 were linked with the G-protein coupled receptor protein signaling pathway, immune response and defense response. DEGs in module 2 were enriched in BPs related to T cell activation and differentiation, and lymphocyte activation and differentiation. DEGs in module 3 were enriched in actin cytoskeleton-related BPs.

Table IV.

GO and KEGG pathway functional enrichment of DEGs in modules.

Category	Term	Description	Count	P-value
Module 1
BP	GO: 0007186	G-protein coupled receptor protein signaling pathway	9	2.20E-09
BP	GO: 0006935	Chemotaxis	6	1.18E-08
BP	GO: 0042330	Taxis	6	1.18E-08
BP	GO: 0007610	Behavior	7	4.44E-08
BP	GO: 0007166	Cell surface receptor linked signal transduction	9	1.24E-07
BP	GO: 0007626	Locomotory behavior	6	1.75E-07
BP	GO: 0006952	Defense response	7	2.23E-07
BP	GO: 0006954	Inflammatory response	6	4.10E-07
BP	GO: 0009611	Response to wounding	6	4.60E-06
BP	GO: 0006955	Immune response	6	1.68E-05
KEGG	hsa04062	Chemokine signaling pathway	8	8.15E-11
KEGG	hsa04060	Cytokine-cytokine receptor interaction	7	1.19E-07
Module 2
BP	GO: 0042110	T cell activation	4	7.78E-06
BP	GO: 0046649	Lymphocyte activation	4	3.07E-05
BP	GO: 0045321	Leukocyte activation	4	5.51E-05
BP	GO: 0001775	Cell activation	4	9.16E-05
BP	GO: 0030217	T cell differentiation	3	2.25E-04
BP	GO: 0030098	Lymphocyte differentiation	3	5.66E-04
BP	GO: 0002521	Leukocyte differentiation	3	9.13E-04
BP	GO: 0007166	Cell surface receptor linked signal transduction	5	1.57E-03
BP	GO: 0045059	Positive thymic T cell selection	2	2.58E-03
BP	GO: 0030097	Hemopoiesis	3	2.93E-03
KEGG	hsa04660	T cell receptor signaling pathway	5	9.47E-07
KEGG	hsa05340	Primary immunodeficiency	3	4.54E-04
KEGG	hsa04640	Hematopoietic cell lineage	3	2.74E-03
KEGG	hsa04514	Cell adhesion molecules	3	6.35E-03
Module 3
BP	GO: 0030036	Actin cytoskeleton organization	4	4.49E-05
BP	GO: 0030029	Actin filament-based process	4	5.44E-05
BP	GO: 0030833	Regulation of actin filament polymerization	3	1.55E-04
BP	GO: 0008064	Regulation of actin polymerization or depolymerization	3	1.98E-04
BP	GO: 0030832	Regulation of actin filament length	3	2.12E-04
BP	GO: 0032271	Regulation of protein polymerization	3	2.47E-04
BP	GO: 0007010	Cytoskeleton organization	4	3.17E-04
BP	GO: 0006928	Cell motion	4	4.08E-04
BP	GO: 0032956	Regulation of actin cytoskeleton organization	3	4.23E-04
BP	GO: 0043254	Regulation of protein complex assembly	3	4.32E-04

GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; DEGs, differentially expressed genes; BP, biological process; Description, the enriched GO function of DEGs; Count, number of DEGs.

Additionally, KEGG pathway analysis for DEGs (Table IV) in module 1 demonstrated that these DEGs predominantly participated in the pathways of chemokine signaling and cytokine-cytokine receptor interaction. DEGs in module 2 participated in T cell receptor signaling pathways, primary immunodeficiency, hematopoietic cell lineage and cell adhesion molecules. However, no pathway was observed to be enriched by the DEGs in module 3.

Data validation demonstrated that the hub genes obtained in GSE15823, such as PTPRC, FN1 and CXCL9, were downregulated. The downregulation of these genes was also observed in GSE41649. Additionally, the expression of CCL21 (enriched in BPs associated with immune response) and NOS2 (enriched in BPs related to oxidation reduction and nitric oxide metabolism) demonstrated the same variation directions in the two datasets.

In addition, following functional enrichment analysis, the DEGs identified in GSE41649 were also demonstrated to be significantly enriched in pathways associated with immune responses (P<0.05), such as systemic lupus erythematosus, complement and coagulation cascades, cytokine-cytokine receptor interaction and chemokine signaling pathways. Furthermore, they were significantly enriched in BPs associated with immune responses and oxidoreductase activity (P<0.05; data not shown).

Discussion

In the present study, the GSE15823 gene expression profile was analyzed and the underlying molecular mechanisms of asthma were explored using bioinformatics methods. A total of 43 upregulated and 275 downregulated DEGs were identified in asthma samples compared with healthy subjects. Upregulated DEGs, such as NOS2, were enriched in BPs related to oxidation, reduction and nitric oxide metabolism. Downregulated DEGs, such as CCL21, CXCL9 and PTPRC, were enriched in immune response-associated BPs, and immune and inflammation-associated pathways. Additionally, 20 DEGs were considered as hub genes in the PPI network, such as PTPRC, CCL21 and CXCL9. The DEGs in module 1 were significantly involved in chemokine signaling pathways. Notably, the 4 DEGs (NOS2, PTPRC, CCL21 and CXCL9), pathways and BP terms were confirmed in GSE41649.

In asthmatic patients, impaired immune response to viral infections is a proposed mechanism for susceptibility to infection (25); however, the relationships between viral infections, innate immunity and asthma are complex and remain to be fully elucidated. In the present study, the majority of the downregulated DEGs were enriched in immune response-associated BPs and pathways, such as CCL21 and CXCL9. Additionally, PTPRC was enriched in BPs related to immune effector processes.

PTPRC was identified as a hub gene with the second highest degree in the PPI network. PTPRC encodes CD45 leukocyte common antigen, which is essential for normal lymphocyte function (26). A study by Trowbridge et al (27) suggested that CD45 has an important role in intracellular signal transduction in the immune system. As such, it has also been demonstrated that the absence of CD45 may result in severe combined immunodeficiency (28). CD45 may be expressed by all leucocytes, including eosinophils (29), which have been implicated as effector cells in asthma and other allergic diseases (30). In asthma, eosinophils and their secreted mediators are believed to be a major contributor to the inflammation underlying the pathophysiological changes in the airways (31). Therefore, we speculated that CD45 (PTPRC) may have an important role in the progression of asthma.

CCL21 and CXCL9 were also identified as hub genes in the PPI network. These genes were also involved in module 1, which was associated with the chemokine signaling pathway. Chemokines are significant regulators of immune cell trafficking, and malfunction of the chemokine signaling pathway has a key role in allergic asthma (32). In allergic asthma, mature dendritic cells migrate to lymph nodes through the interaction of chemokine (C-C motif) receptor 7 and CC chemokine ligands of CCL21, where they effectively present antigens to naive T cells (33). It has previously been demonstrated that the administration of antigen-pulsed dendritic cells induces allergic inflammation (34). CCL21 is understood to be crucial for lymphoid cell trafficking and for the structural organization of lymphoid tissues (35). A study by Xu et al (36) suggested that a lack of lymphoid chemokine ligand, CCL21, enhances allergic airway inflammation by modulating the recruitment of CD4+ T cells into the lungs. CXCL9 may be produced by eosinophils following stimulation with interferon-γ, indicating that CXCL9 may be involved in the downregulation of allergic inflammation (37,38). Taken together these findings suggest that CCL21 and CXCL9 may be involved in immune response-related BPs and have important roles in the progression of asthma.

In addition to immune responses, oxidants may also have a critical role in the pathogenesis of asthma (39). The lungs are highly susceptible to oxidative injury, which has been implicated in various lung diseases, including asthma (40). Increased levels of reactive oxygen species may induce apoptosis and result in increased airway reactivity and secretions, which may augment existing inflammation in asthma (41). In the present study, BPs related to oxidation reduction and nitric oxide metabolism were enriched by several upregulated DEGs, including NOS2. NOS2 is able to endogenously produce NO, which has been suggested to have an important role in the physiological regulation of airway functions (42,43). Notably, a previous study demonstrated that NOS2 was significantly associated with the fraction of exhaled NO in asthmatic children (44). Additionally, a study by Dweik (45) reported an increase in NOS2 expression during inflammation, which was in accordance with the results of the present study. Taken together these findings suggest that NOS2 may be a candidate molecular marker associated with asthma progression through oxidation reduction.

In the present study, there were 4 lung tissue samples obtained from bronchial biopsies of healthy controls and 4 from subjects with allergic asthma. All 4 of the asthmatic patients were women, whereas the 4 healthy subjects included 2 women and 2 men. Notably, Laprise et al (10) stated that gene expression measurements were similar in control subjects, which indicated that there was no gender difference in our data. Similarly, a study by Prescott et al (46) demonstrated that there was no gender difference in the prevalence of asthma.

The dataset of GSE41649 confirmed the key genes, pathways and GO functions identified in the present study; however, some limitations existed. In the process of data analysis, only 8 samples were used. Furthermore, the investigation was implemented by means of bioinformatics methods and the screened genes and pathways have not been validated by experiments, although the dataset of GSE41649 was used for verification.

In conclusion, the results of the present study demonstrate that the immune response, oxidants and nitric oxide metabolism may have important roles in the progression of asthma. Additionally, DEGs, such as PTPRC, CXCL9, CCL21 and NOS2, may have the potential to be used as targets for asthma diagnosis and treatment.