Protein–protein interaction network analysis and identifying regulation microRNAs in asthmatic children

Abstract

Background

Asthma is a chronic inflammatory airway disease, the incidence of which has increased recently. In order to identify the potential biomarkers for allergic asthma therapy, microarray data were analysed to find meaningful information.

Methods

Microarray data GSE18965 were downloaded from Gene Expression Ominibus (GEO), including seven asthmatic epithelium samples from children with allergic asthma and nine healthy controls. Limma package was used to detect differentially expressed genes (DEGs) and the criteria were |log fold change| > 0.5 and p value < 0.05. We used Database for Annotation, Visualization and Integrated Discovery (DAVID) tool to perform GO function and KEGG pathway analysis. STRING database was used to construct protein–protein interaction (PPI) network. MicroRNA (miRNA) regulation network was constructed according to miRecords database.

Results

We identified 274 DEGs in asthma epithelium samples comparing with healthy controls. There were 123 up-regulated DEGs and 151 down-regulated DEGs. PPI network analysis showed that TSPO, G6PD and TXN had higher degree. miRNA regulation network demonstrated that miR-16 and miR-15a had higher degree. The target genes of miRNAs were significantly enriched in the apoptosis function.

Conclusions

TSPO, G6PD and TXN, miR-16, miR-15a and apoptosis may be used as the targets for children's allergic asthma therapy.

Introduction

Asthma is a common chronic inflammatory airway disease, which is characterised by variable airflow obstruction, hyperresponsiveness, airway inflammation and reversibility either spontaneously or as a result of treatment.¹ Asthma is triggered by virus and pollutants through exposure of airway epithelial cell (AEC). The damage of AEC results in increased epithelial permeability and penetration of allergens and injurious materials, reflecting an inadequate repair response, an increased susceptibility to injury, or a combination of both.² There have been a total of 334 million individuals suffering from asthma worldwide. Asthma is the most common chronic disease for children, and death rates range from 0.0 to 0.7 per 100 000 in children.^3,4 Thus, it is necessary to explore the mechanism of asthma to achieve effective and satisfactory control of asthma in children.

Asthma is characterised by airway hyperreactivity, which is mediated by Th2 lymphocytes producing interleukin-4 (IL-4), IL-5, IL-9 and IL-13.⁵ T cell- and eosinophil-induced apoptosis has been reported to be an important event in asthma.⁶ In severe asthma, there is increased cellular proliferation in the airway contributing to a thickened epithelium and lamina reticularis.⁷ Microarray profiling has been widely used for the analysis of expressed genes associated with asthma. Woodruff et al. found that calcium-activated Cl-channel-1 and plasminogen activator inhibitor-2 were up-regulated in asthma by genome-wide profiling.⁸ Fibronectin was the only extracellular matrix component, the expression of which was significantly lower in asthmatic AEC.⁹

Several reports have reported that microRNAs (miRNAs) played important roles in the pathology of asthma. miRNA is small non-coding RNAs and plays a regulatory role.³ miRNA-221 was differentially expressed in asthmatics and controls, suppression of which could inhibit airway inflammation in asthmatics.¹⁰ Toll-like receptor 7 was significantly reduced in severe asthma compared to human controls, which was regulated by miRNAs.¹¹ miRNA-221 and miRNA-485-3p have been reported to be up-regulated in paediatric asthma.¹²

In this study, differentially expressed genes (DEGs) in the allergic asthma samples comparing with healthy controls were screened out. Then Gene Ontology (GO) function and kyoto encyclopaedia of genes and genomes (KEGG) pathway analysis were performed for the DEGs. Subsequently, we constructed protein–protein interaction (PPI) network of these DEGs. Besides, regulating network of interactions between miRNAs and DEGS was constructed according to the information in the miRecords dataset.

Materials and methods

Microarray data

Gene expression profiling GSE18965 of airway epithelial cells was downloaded from Gene Expression Ominibus (GEO) including seven allergic asthma samples and nine healthy controls from children. The annotation platform was GPL96 [HG-U133A] Affymetrix Human Genome U133A Array.

Affymetrix microarray CEL files were converted into expression profile matrix by Robust Multi-array Average (RMA) algorithm of affy package in the R language.¹³ Probes were mapped to gene name according to annotation information in the platform Affymetrix by Bioconductor software.¹⁴ Expression values of multiple probes that correspond to one gene were averaged.

DEGs screening and function analysis

DEGs were screened by linear models for microarray data (limma) package in R language.¹⁵ The criteria were |log FC| > 0.5 (FC: fold change) and p < 0.05. Database for Annotation, Visualization and Integrated Discovery (DAVID) online tool was used for GO and KEGG pathway analysis with p < 0.05.^16–18

PPI network analysis

STRING online database was used to construct PPI network of DEGs.¹⁹ The hub genes were screened according to the degree of node. The degree represents the number of interaction partners and was calculated by Perl code.²⁰ Then the cytoscape was used to construct sub-network.²¹ Then DAVID tool was used to perform function analysis for the sub-network module.

Regulation relation between miRNA and DEGs

The information of miRNA was extracted from miRecords²² and then the miRNA regulation network was constructed. The interactions of miRNAs and targets genes have been experimentally validated previously. Finally, target genes of miRNAs performed GO function analysis.

Results

DEGs in allergic asthma samples comparing with healthy controls

A total of 12 495 DEGs were obtained after normalisation. Then, 274 DEGs were identified in the asthma samples comparing with healthy controls. Among them, 123 DEGs were up-regulated and 151 genes were down-regulated.

GO function and KEGG pathway analysis

We categorised DEGs by known or proposed gene function (Table 1 listed the top 19 categories). These up-regulated genes participated in the regulation of transcription, negative regulation of nucleic acid and nitrogen compound metabolic process; were involved in the cellular components that are intracellular non-membrane-bounded organelle, non-membrane-bounded organelle, and cytosol; and were involved in molecular functions that are RNA combination, structure of ribosome and myosin light streptokinase activity, as well as pathways that may be specifically related to regulation of actin cytoskeleton regulation and ribosomes. Down-regulated genes participated in the response of the organic substance, protein localisation and intracellular transport; were involved in the cellular components such as plasma membrane and cytosol; and were involved in molecular functions such as peptidase activity, endopeptidase activity and actin binding activity, as well as pathways that may be specifically related to metabolism of lysosomes, glutathione and galactose.

Table 1.

GO and KEGG enrichment analysis for differentially expressed genes (DEGs). The top 10 BP (biological process), MF (molecular function), cellular component (CC) terms with DEGs which were significantly enriched were listed.

Genes	Category	Term	Count	p value
Up-regulate genes	GOTERM_BP_FAT	GO: 0045449∼regulation of transcription	28	0.021779
	GOTERM_BP_FAT	GO: 0006350∼transcription	24	0.019963
	GOTERM_BP_FAT	GO: 0045934∼negative regulation of nucleobase, nucleoside, nucleotide and nucleic acid metabolic process	9	0.028924
	GOTERM_BP_FAT	GO: 0051172∼negative regulation of nitrogen compound metabolic process	9	0.031005
	GOTERM_BP_FAT	GO: 0051253∼negative regulation of RNA metabolic process	8	0.014452
	GOTERM_BP_FAT	GO: 0016481∼negative regulation of transcription	8	0.044573
	GOTERM_BP_FAT	GO: 0043009∼chordate embryonic development	7	0.030307
	GOTERM_BP_FAT	GO: 0009792∼embryonic development ending in birth or egg hatching	7	0.031476
	GOTERM_BP_FAT	GO: 0045892∼negative regulation of transcription, DNA-dependent	7	0.04096
	GOTERM_BP_FAT	GO: 0030036∼actin cytoskeleton organisation	6	0.02232
	GOTERM_CC_FAT	GO: 0043232∼intracellular non-membrane-bounded organelle	26	0.009712
	GOTERM_CC_FAT	GO: 0043228∼non-membrane-bounded organelle	26	0.009712
	GOTERM_CC_FAT	GO: 0005829∼cytosol	16	0.01356
	GOTERM_CC_FAT	GO: 0005794∼Golgi apparatus	12	0.016346
	GOTERM_CC_FAT	GO: 0030529∼ribonucleoprotein complex	8	0.037376
	GOTERM_CC_FAT	GO: 0044445∼cytosolic part	6	0.00234
	GOTERM_CC_FAT	GO: 0005840∼ribosome	6	0.010051
	GOTERM_CC_FAT	GO: 0033279∼ribosomal subunit	5	0.007734
	GOTERM_CC_FAT	GO: 0015934∼large ribosomal subunit	4	0.007914
	GOTERM_CC_FAT	GO: 0022626∼cytosolic ribosome	4	0.013265
	GOTERM_MF_FAT	GO: 0003723∼RNA binding	11	0.028442
	GOTERM_MF_FAT	GO: 0003735∼structural constituent of ribosome	5	0.030331
	GOTERM_MF_FAT	GO: 0004687∼myosin light chain kinase activity	2	0.020882
	KEGG_PATHWAY	hsa04810: regulation of actin cytoskeleton	6	0.027998
	KEGG_PATHWAY	hsa03010: ribosome	4	0.032319
Down-regulate gene	GOTERM_BP_FAT	GO: 0010033∼response to organic substance	17	0.00137
	GOTERM_BP_FAT	GO: 0008104∼protein localisation	17	0.009684
	GOTERM_BP_FAT	GO: 0046907∼intracellular transport	15	0.00396
	GOTERM_BP_FAT	GO: 0055114∼oxidation reduction	14	0.007952
	GOTERM_BP_FAT	GO: 0006955∼immune response	14	0.014456
	GOTERM_BP_FAT	GO: 0042592∼homeostatic process	14	0.026908
	GOTERM_BP_FAT	GO: 0045184∼establishment of protein localisation	14	0.031774
	GOTERM_BP_FAT	GO: 0006091∼generation of precursor metabolites and energy	12	2.01E-04
	GOTERM_BP_FAT	GO: 0009719∼response to endogenous stimulus	11	0.005385
	GOTERM_BP_FAT	GO: 0009725∼response to hormone stimulus	10	0.008521
	GOTERM_CC_FAT	GO: 0005886∼plasma membrane	53	0.017878
	GOTERM_CC_FAT	GO: 0044459∼plasma membrane part	35	0.013535
	GOTERM_CC_FAT	GO: 0005829∼cytosol	26	0.003322
	GOTERM_CC_FAT	GO: 0000267∼cell fraction	22	0.0051
	GOTERM_CC_FAT	GO: 0005739∼mitochondrion	22	0.005318
	GOTERM_CC_FAT	GO: 0031226∼intrinsic to plasma membrane	21	0.032332
	GOTERM_CC_FAT	GO: 0005887∼integral to plasma membrane	20	0.046866
	GOTERM_CC_FAT	GO: 0005783∼endoplasmic reticulum	19	0.013103
	GOTERM_CC_FAT	GO: 0005773∼vacuole	18	1.37E−09
	GOTERM_CC_FAT	GO: 0005794∼Golgi apparatus	18	0.011081
	GOTERM_MF_FAT	GO: 0070011∼peptidase activity, acting on l-amino acid peptides	12	0.011274
	GOTERM_MF_FAT	GO: 0008233∼peptidase activity	12	0.015298
	GOTERM_MF_FAT	GO: 0004175∼endopeptidase activity	9	0.021001
	GOTERM_MF_FAT	GO: 0003779∼actin binding	8	0.029168
	GOTERM_MF_FAT	GO: 0048037∼cofactor binding	7	0.026057
	GOTERM_MF_FAT	GO: 0050662∼coenzyme binding	6	0.024665
	GOTERM_MF_FAT	GO: 0003924∼GTPase activity	6	0.043348
	GOTERM_MF_FAT	GO: 0050661∼NADP or NADPH binding	4	0.0042
	GOTERM_MF_FAT	GO: 0051015∼actin filament binding	4	0.012327
	GOTERM_MF_FAT	GO: 0015078∼hydrogen ion transmembrane transporter activity	4	0.048639
	KEGG_PATHWAY	hsa04142: lysosome	13	4.32E−08
	KEGG_PATHWAY	hsa00480: glutathione metabolism	6	5.61E−04
	KEGG_PATHWAY	hsa00052: galactose metabolism	5	3.92E−04
	KEGG_PATHWAY	hsa00030: pentose phosphate pathway	4	0.004625
	KEGG_PATHWAY	hsa00520: amino sugar and nucleotide sugar metabolism	4	0.022122
	KEGG_PATHWAY	hsa00010: glycolysis/gluconeogenesis	4	0.048967
	KEGG_PATHWAY	hsa00410: beta-alanine metabolism	3	0.03614

PPI network analysis

Fig. 1 showed the PPI network of DEGs, demonstrating that TSPO (translocator protein), G6PD (glucose-6-phosphate dehydrogenase deficiency) and TXN (thioredoxin) have the higher degrees which were 38, 20 and 17, respectively. The interactions among these three nodes and their first adjacent nodes were used to construct the sub-network (Fig. 2). The major biology processes which the sub-network modules participated in were homeostatic process, protein localisation and regulation of cell proliferation (Table 2).

Figure 1.

Protein–protein interaction network which the differentially expressed genes corresponded to. Red dots represent up-regulated genes and green dots represent down-regulated genes. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of the article.)

Figure 2.

Sub-network of top three hub proteins and their first adjacent nodes of hub proteins. Red dots represent up-regulated genes and green dots represent down-regulated genes. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of the article.)

Table 2.

Top 10 biology processes in which genes in the sub-network were significantly enriched.

Term	Count	p value
GO: 0042592∼homeostatic process	15	9.53E−07
GO: 0008104∼protein localisation	13	1.47E−04
GO: 0042127∼regulation of cell proliferation	12	2.31E−04
GO: 0019725∼cellular homeostasis	11	1.31E−05
GO: 0055114∼oxidation reduction	10	8.74E−04
GO: 0046907∼intracellular transport	10	0.001063
GO: 0008219∼cell death	10	0.00199
GO: 0016265∼death	10	0.002086
GO: 0045184∼establishment of protein localisation	10	0.003136
GO: 0006091∼generation of precursor metabolites and energy	9	3.09E−05

Regulation network of miRNA and DEGs

Fig. 3 shows the relation network between miRNA and DEGs. In the network, miR-16 and miR-15a have the higher degrees, which were 62 and 59, respectively. The target genes of miRNA participated in the cell apoptosis, programmed cell death and cell death (Table 3).

Figure 3.

Interaction network of miRNAs and differentially expressed genes (DEGs). Yellow dots represent miRNAs and green dots represent DEGs. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of the article.)

Table 3.

Top 10 biology processes in which the targets genes of miRNA were enriched.

Term	Count	p value
GO: 0006915∼apoptosis	14	3.14E−05
GO: 0012501∼programmed cell death	14	3.66E−05
GO: 0008219∼cell death	14	1.90E−04
GO: 0016265∼death	14	2.04E−04
GO: 0043933∼macromolecular complex subunit organisation	13	6.32E−04
GO: 0033554∼cellular response to stress	12	3.39E−04
GO: 0065003∼macromolecular complex assembly	12	0.001291
GO: 0006461∼protein complex assembly	11	5.58E−04
GO: 0070271∼protein complex biogenesis	11	5.58E−04
GO: 0006259∼DNA metabolic process	10	0.002246

Discussion

In this study, we identified 274 DEGs in the asthma samples comparing with healthy controls. Then we used PPI network to identify hub genes and screened key miRNAs which may be associated with asthma.

PPI network analysis showed that TSPO, G6PD and TXN had the high degree. TSPO is an 18 kDa protein mainly found in the mitochondrial membrane,²³ which has been shown to alter mitochondrial respiration in a dose dependent manner, increasing the rate of succinate-supported state IV respiration and decreasing ADP-supported state III respiration, thereby decreasing the respiratory control ratio.²⁴ The severity of atopic dermatitis, a chronic inflammatory skin disease, was related with TSPO expression in male patients.²⁵ It has been reported that TSPO is highly expressed on the bronchial and bronchiole epithelium in the lungs of humans and guinea pigs,²⁶ pneumocytes and alveolar macrophages²⁷ in humans. However, there were no reports about the relation between TSPO and asthma. In this study, TSPO had a high degree in the PPI network. The hub protein which had a high degree tends to be essential.²⁸ Therefore, we speculated that TSPO might be a potential biomarker for asthma therapy and provided basic understanding for the pathologic understanding of asthma. G6PD is a cytosolic enzyme that catalyses the chemical reaction from d-glucose 6-phosphate to 6-phospho-d-glucono-1,5-lactone.²⁹ In individual embryos, mRNAs expression of G6PD was correlated with respiration rates and uptake of oxygen.³⁰ The expression of antioxidant gene G6PD was induced in the enhanced severity of the asthmatic mice in response to allergen through RT-PCR analysis.³¹ In acute exacerbations of asthma, the serum TRX levels in patients with asthma were significantly increased compared to those during the asymptomatic period.³² There was a high sensitisation rate of thioredoxin among bakers.³³ In ovalbumin-sensitised mice, TRX strongly suppressed airway hyperresponsiveness and airway inflammation.³⁴ Thus, G6PD and TRX may be biomarkers in the airway epithelial cells for asthma therapy.

The sub-network of hub proteins (TSPO, G6PD and TRX) was significantly enriched in the homeostatic process. Surfactant homeostasis was altered in the Aspergillus fumigatus-induced allergic airway inflammation, demonstrating the down-regulation of hydrophobic surfactant proteins (SPs), SP-B and SP-C and increase in SP-D protein levels.³⁵ Altered l-arginine homeostasis and deficiency of nitric oxide due to increased arginase activity have been considered to be the causative factors in the pathology of asthma.³⁶ Thus, the homeostasis process may be one of the important biological processes in the pathogenesis of allergic asthma.

Interaction network of miRNA and DEGs showed that miR-16 and miR-15a had high degree. Lower levels of miR-15a decreased the expression of vascular endothelial growth factor-A which is over-expressed in the sputum and serum of both adult and paediatric patients with asthma.³⁷ miR-16 was up-regulated in the mice when exposing the airways to house dust mite by quantitative PCR analysis.³⁸ The miR-15 and miR-16 in the network demonstrated their relation with the DEGs in the asthma tissue. We could conclude that miR-15 and miR-16 might regulate the genes in the asthma tissue in response to stimulus. The target genes of miRNAs were significantly enriched in apoptosis, programmed cell death, and cell death. Asthma is characterised by the accumulation of eosinophils into the lungs, which participate in the maintenance and exacerbation of inflammation. Inflamed airway eosinophil apoptosis has been observed by anti-asthmatic agent such as glucocorticoids, theophylline and leukotriene modifiers, which may be a therapeutic target in the treatment of allergic asthma.³⁹ Besides, it has been reported that the spontaneous and mitogen-induced apoptosis of mononuclear cells were increased in patients with atopic asthma.⁴⁰ Thus, apoptosis, also called programmed cell death, may be the potential target for allergic asthma.

In summary, TSPO, G6PD and TXN with a high degree in the PPI network may be potential targets for asthma therapy. In addition, miR-15 and miR-16 might be the important miRNAs which regulate the DEGs in allergic asthma. The apoptosis pathway could also be the target for asthma therapy. Lack of experiments is the drawback of our study. So the experiments for investigating the roles of TSPO, G6PD, TXN, miR-15 and miR-16 in the allergic asthma would be needed in a future study.

Conflict of interest

The authors declare that they have no conflict of interests to state.

Authors’ contributions

Y.Q.C. carried out the design and coordinated the study, participated in most of the experiments and prepared the manuscript. J.N.Q. provided assistance in the design of the study, coordinated and carried out all the experiments and participated in manuscript preparation. All authors have read and approved the content of the manuscript.

Funding

The National Natural Science Fund Project (No. 81170028).

Ethical disclosures

Patients’ data protection

Confidentiality of data. The authors declare that they have followed the protocols of their work centre on the publication of patient data and that all the patients included in the study have received sufficient information and have given their informed consent in writing to participate in that study.

Right to privacy and informed consent

The authors have obtained the informed consent of the patients and/or subjects mentioned in the article. The author for correspondence is in possession of this document.

Protection of human and animal subjects

The authors declare that the procedures followed were in accordance with the regulations of the responsible Clinical Research Ethics Committee and in accordance with those of the World Medical Association and the Helsinki Declaration.