Transcriptome analysis of peripheral whole blood identifies crucial lncRNAs implicated in childhood asthma

Background Asthma is a chronic disorder of both adults and children affecting more than 300 million people heath worldwide. Diagnose and treatment for asthma, particularly in childhood asthma have always remained a great challenge because of its complex pathogenesis and multiple triggers, such as allergen, viral infection, tobacco smoke, dust, etc. It is thereby great significant to deeply investigate the transcriptome changes in asthmatic children before and after desensitization treatment, in order that we could identify potential and key mRNAs and lncRNAs which might be considered as useful RNA molecules for observing and supervising desensitization therapy for asthma, which might guide the diagnose and therapy in childhood asthma. Methods In the present study, we performed a systematic transcriptome analysis based on the deep RNA sequencing of ten asthmatic children before and after desensitization treatment, including identification of lncRNAs using a stringent filtering pipeline, differential expression analysis and network analysis, etc. Results First, a large number of lncRNAs were identified and characterized. Then differential expression analysis revealed 39 mRNAs and 15 lncRNAs significantly differentially expressed which involved in two biological processes and pathways. A co-expressed network analysis figured out a desensitization-treatment-related module which contains 27 mRNAs and 21 lncRNAs using WGCNA R package. Module analysis disclosed 17 genes associated to asthma at distinct level. Subsequent network analysis based on PCC figured out several key lncRNAs probably interacted to those key asthma-related genes, i.e., LINC02145, GUSBP2. Our functional investigation indicated that their functions might involve in immune, inflammatory response and apoptosis process. Conclusions Our study successfully discovered many key noncoding RNA molecules related to pathogenesis of asthma and relevant treatment, which may provide some clues for asthmatic diagnose and therapy in future.

prevalence rates of asthma among children increases from 5% in Albania to 21% in the UK. 9.1% of US children (6,7 million) had asthma in 2007 [1]. In China, the prevalence of childhood asthma increased extremely since 1990s, which ranged from 0.93% in 1990 to 1.54% in 2000 [3]. There is no definitive explanation for this increase, even though the prevalence of asthma in children has increased over decades. More importantly, there is no full cure is available up to now despite majority of children with asthma can alleviate the asthmatic symptoms and obtain adequate asthma control via avoidance of triggering factors, rational management and/or medication, such as short-acting inhaled β2-receptor agonists [4,5]. However, a small proportion (~5%) of asthmatic children have uncontrolled asthma despite maximum medical treatment [6]. In addition, diagnosing asthma in children faces challenge as well, e.g., a number of childhood conditions exhibit compatible symptoms to those caused by asthma, such as shortness of breath, wheezing and cough [5]. To further complicate the issue, those conditions can coexist with asthma and confuse the evaluation of patients.
Like other immune-related diseases, genetic components, i.e., sequences variants, in combination with environmental factors, i.e., allergen, viral infection, may contribute to the development of asthma [7][8][9]. Concretely, genome-wide association studies have detected numerous asthma-associated gene variants, while few of them (less than 10%) can be explained to contribute to the risk of asthma [10]. On the other hand, environmental factors have been reported to modulate epigenetic modifications and thereby affect gene expression and phenotype [11,12]. Epigenetic mechanisms have been reported to be key to functions of airway smooth muscle (ASM), e.g., aberrant phenotype of ASM lead to the obstruction and inflammation of airway [13]. Therefore, diverse cellular compartments orchestrated by multiple environment-driving factors make the mechanisms underlying asthma extraordinary complex and unclear, particularly with respect to long noncoding RNAs (lncRNAs).
lncRNAs have been reported to possess multiple functions including regulation of gene expression, transcriptional activation and silencing of genes and thus play critical roles in a diversity of cellular process [14,15]. Aberrant expression of many specific lncRNAs was reported to be correlated to the pathogenesis of various diseases in human and model animals [16,17]. Furthermore, the roles of lncRNAs implicated in asthma have been emerged recently. Several studies achieved by whole genome RNA-seq have identified many lncRNAs specifically expressed in T cells at varying stages of development and differentiation [18]. For instances, the lncRNA named Tmevpg1 was found to be uniquely expressed by the T helper cell 1 (TH1) and is essential to the transcription of Ifng by TH1 subset [19]. The lncRNA GATA3-AS1 transcribed in the opposite direction from GATA3 was reported to be co-expressed with GATA3 in mouse and human TH2 cells [20]. TH2 cells responses have already demonstrated to link to both allergy and asthma. However, the gaps between the function of the lncRNAs specifically expressed in immune system or T lymphocytes and asthma still unfilled and require further deep investigation, e.g., despite GATA3-AS1 may play a role in asthmatic response, its functions are unknown. Here, to attempt to fill some of those gaps and identify the novel lncRNAs which might be useful molecules for observing and supervising desensitization therapy for childhood asthma, we performed a systematic transcriptome analysis based on the deep RNA-seq data of peripheral whole blood of ten asthmatic children before and after desensitization treatment.

Identification and characteristics of lncRNAs
At first, the raw data was filtered and trimmed to eliminate the adaptor sequences and low-quality reads using Trimmomatic. This process yielded an average of 121, 690,584 clean reads for each sample (Table S1). Then the clean reads were utilized to re-construct the new transcriptome for human for the present study. Firstly, all the clean data were aligned to human genome sequences using STAR (Table S2). After that, the aligned bam files were subjected to StringTie to assemble into 215,331transcripts (Table S3). We simply compared these assembled transcripts against the known transcripts of human derived from Ensembl database via BLASTn, the result showed that 16,220 (~7.5%) novel transcripts were identified in the present study.
To identify high-confidence lncRNAs for this study, we initially extracted the known lncRNAs based on the annotation information of gtf file of Homo sapiens (GRCh38.p10) from Ensembl database (details see Methods). As a result, a total of 52,791 known lncRNAs were detected. On the other hand, we applied a stringent stepwise filtering pipeline ( Fig. 1) to identify novel lncRNAs. This analysis generated 5310 novel lncRNAs (Table S4). Furthermore, aimed at assessing the quality of identified lncRNAs, we performed a comparative study between lncRNAs and mRNAs before and after treatment in the transcript length, exon counts, open reading frame (ORF) length and expression level (Fig. 2). The results indicated that there was no obvious difference of expression level among these patients no matter in mRNAs or lncRNAs (Fig. 2a). The similar situation exists in the comparison of expression before and after treatment (Fig. 2b). In consistent patterns of with previous studies of eukaryotic lncRNAs, the length of the majority of identified lncRNAs and the corresponding ORF was much shorter than that of mRNAs (Fig. 2c, d), and the exon count of the lncRNAs was also less than that of mRNAs (Fig. 2e).
Moreover, we found that a small proportion of known and novel lncRNAs specifically expressed in distinct status (Fig. 3a, b). In addition, the length of majority of known lncRNAs ranges from 500 bp to 700 bp, whereas the length of lncRNAs were found mainly shorter than 300 bp (Fig. 3c). We also detected five main alternative splicing (AS) modes for these transcripts using AStalavista [21]. The results showed that there was no significant discrepancy in those five AS modes before and after treatment (Fig. 3d).

Differential expression analysis reveals desensitizationtreatment-related module
To investigate the lncRNAs as well as their possible roles implicated in asthma before and after treatment, we initially performed paired differential-expression analysis using DESeq2 package for R [22]. The results indicated that only 54 transcripts (15 lncRNAs and 39 mRNAs) were significantly differentially expressed in the patients after treatment (Table 1, Fig. 4a). After that, all the differentially expressed transcripts were subjected to metascape to search for possible functional biological processes. As a result, two biological processes were significantly enriched, including humoral immune response and cytoplasmic vesicle membrane (Fig. 4b).
Co-expression network analysis reveals desensitizationtreatment-related module Alternatively, we utilized a more sensitive and robust approach to explore the lncRNAs implicated in asthma, i.e., co-expression network analysis, which was achieved via WGCNA in R. Firstly, the expression profile of all transcripts was assessed using featurecounts and normalized by DESeq2 package for R. Then the cluster analysis was conducted using flash tools package of WGCNA to detect outliers of samples. The results indicated that all samples were well clustered and passed the cuts ( Figure S1). A network-topology analysis based on different soft- Fig. 1 The bioinformatics pipeline used for the identification of lncRNAs for ten asthmatic children thresholding powers was performed to get relative balance scale independence and mean connectivity of WGCNA. The optimal power for the scale-free topology fit index was calculated out for construction of hierarchical clustering tree ( Figure S2). This analysis clustered all the correlated transcripts into 80 modules with strong correlation with trait ( Figure S3). Each module contains independent datasets of transcripts (includes mRNAs and lncRNAs). The interactions among those modules were visualized in Figure S4 with randomly selected transcripts. After that, the modules significantly associated with specific trait (desensitization treatment) was detected using the function plotEigengeneNetworks of WGCNA ( Table 2). The results showed that a module (khaki1) were significantly correlated to the trait (desensitization treatment) (Fig. 5), which contains 27 mRNAs and 21 lncRNAs. The relationship of transcript significance with module membership was visualized in Figure S5.  2 General characteristics of assembled transcripts in ten asthmatic children between two different conditions. a Expression level indicated by log 10 (reads count) in mRNAs and lncRNAs in for ten asthmatic children before and after treatment. b Expression level indicated by log 10 (reads count) in lncRNAs between two different conditions. c Distribution of transcript length by log 10 (Length) d ORF length by log 10 (ORF Length) e Exon count by log 10 (Exon Count) in mRNAs and lncRNAs in in ten asthmatic children Functional analysis of module reveals key genes related to asthma In order to investigate the potential relationship of transcripts which co-localizes in desensitization-treatmentrelated module with asthma treatment, the annotated transcripts (mRNAs) of the module were subjected to a gene annotation server Metascape for functional enrichment analysis. The result showed that one biological process (GO:0043065: positive regulation of apoptotic process) was significantly enriched ( In addition to these four apoptosis-related genes (TRIO, PEA15, KLRK1, NDUFA13), we also found that the other genes in desensitization-treatment-related module  exhibiting distinct degree of correlations to asthma via literature review [23][24][25][26][27][28][29][30][31][32][33][34][35][36][37]. Concretely, four genes (ATF7IP, SLC43A3, AKAP7, HMGN1) were found to be correlated to pathogenesis or treatment of asthma in immune. For instance, ATF7IP can form a protein-complex by binding to MBD1 (methyl CpG binding protein 1) and Aire (autoimmune regulator) to maintain immune tolerance [23]. Immunological tolerance to self-antigens is critical to the prevention of autoimmune disease, including asthma. In addition, two genes (ELMO2, WDFY1) were demonstrated to play roles in inflammation response of asthma. ELMO2 was found as an ILK-binding protein to play key role in integration of adhesion and migration pathways. Abnormal regulation of migration is associated to chronic autoinflammatory disorder, such as asthma [33]. Furthermore, two genes (G2E3, ZNF19) were found to be differentially methylated between asthmatics and nonasthmatics [34,35]. Combined with four apoptosis-related genes, all those asthma-related genes could be summarized in Table 3.

Network analysis based on PCC reveals key lncRNAs implicated in asthma
Our aforementioned analysis has successfully assisted us to identify several asthma-related genes in desensitizationtreatment-related module. However, the correlation of lncRNAs with asthma in the module were still unknown. For that, we performed a network analysis to establish RNA-RNA interaction network for all the 48 transcripts (including mRNAs and lncRNAs) in desensitizationtreatment-related module based on the remaining 102,041 assembled transcripts using pearson correlation coefficient (PCC) method. The RNA-RNA pairs which had absolute value of pearson correlation ≥0.7 and P-value< 0.05 were retained. This analysis detected a total of 3502 RNA-RNA interaction relationships including mRNA-mRNA, mRNA-lncRNA, which involved 1654 mRNAs and 864 lncRNAs (Fig. 6).
We further concerned which lncRNAs target to the 17 asthma-related genes as well as their possible roles in asthma. Since that, we extracted a sub-network of these 17 asthma-related genes from the above network (Fig. 7a). This step indicated that 323 lncRNAs exhibited potential correlations to these asthma-related genes, implying they may play roles in pathogenesis and treatment of asthma. To further investigate their possible functions implicated in asthma, the target mRNAs for each asthma-related gene were extracted to be subjected to Metascape server for functional analysis. The results showed that majority of target mRNAs of asthma-related gene could significantly enriched by diverse biological processes (Fig. 7a).
Among them, many RNA-RNA interactions aroused our attentions since both sides of interaction co-localizes in the desensitization-treatment-related module (Fig. 7b), which suggested their strong correlation with asthma. Concretely, three known lncRNAs (LINC01771, LINC02145, AL031289.1) and three novel lncRNAs (MSTRG.24219.1, MSTRG.28178.1, MSTRG.18506.1) were found to interact with many asthma-related genes. Notably, LINC02145 exhibited to be correlated to six asthma-related genes, including ZFN19, G2E3, ATF7IP, PEA15, SPDYE6, VWA5A. Despites there is no clinical significance was reported in this known lncRNA, our analysis indicated that its strong association with asthma before and after desensitization treatment. A similar situation was found in a novel lncRNA MSTR G.18506.1 which interacts with three asthma-related genes, i.e., SPDYE6, ATF7IP, VPS37A. Interestingly, in addition to lncRNAs, we found some other types of noncoding RNAs exhibited correlations with asthma-related genes, including antisense, snoRNA, snRNA and pseudogene. To test our findings that those specific lncRNAs and mRNAs probably associated to immune/inflammatory-related genes in asthma, we selected three lncRNA/mRNAs, including KLRK1, LINC02145, GUSBP2, which exhibits differential expression before and after treatment for the qRT-PCR. The results were demonstrated to be consistent to our bioinformatic analysis ( Figure S7).

Discussion
With the advancement of high throughput sequencing technology, tens of thousands of lncRNAs have been detected in human. Despite the function of majority of Fig. 4 Differential expression analysis between before treatment and after treatment. a Heatmap of expression profiles using differentially expressed transcripts. The heatmap was visualized using R package gplots, all expression values were normalized using Z-score using R package DESeq2. b Functional GO enrichment analysis based on differentially expressed transcripts using Metascape lncRNAs are unknown so far, their important roles in diverse biological processes has already been demonstrated, particularly close correlations to diseases. In the present study, we detected at least one order of magnitude lncRNAs larger than those identified in our previous studies about lncRNAs [38,42,43]. One of reason is the strategy of RNA sequencing in our previous studies is to capture the RNAs with polyA tail (mRNAs), which could only identify around half of lncRNAs. Here we applied total RNA sequencing technology to capture all linear RNAs, which contains most of lncRNAs. In addition, we improved the analysis pipeline for identification of lncRNAs used in our previous studies, which mainly focused on the identification of novel lncRNAs. Then the differential expression analysis showed that only a small proportion of RNAs were differentially expressed in asthmatic patients before and after treatment. The possible reason is that asthma as a chronic respiratory disease, treatment of it might not directly be reflected in peripheral blood. The RNA change triggered by treatment in asthmatic patients might mainly occur in airway epithelium. This is one limitation of our study. This limitation could be solved by using the deep RNA sequencing data of airway epithelium cell. Another limitation to be noted is that ten individuals participated in the present study are composed of nine male children and one female child. To our knowledge, gender, in particular during pubertal period (one female child in this study is 15 years old) might affect the results of differential expression analysis. To be clear, the results and the conclusions in the present study are based mainly on male population of asthmatic childhood patients. Regardless, the functional enrichment analysis of differentially expressed transcripts indicated their functions might be involved in humoral immune response. Humoral immune response is one of two main arms of the immune system, which triggers specific B cells to proliferate and secrete large amounts of their specific antibodies. This process attracts a helper T (T H ) cell to be activated against a particular antigen, which implying its importance in immune response to asthma.
After that, co-expression network analysis was achieved by WGCNA package to identify desensitizationtreatment-related module. Although correlation value of the most trait-related (khaki1) module is not very significant (only reach 0.63) compared to other WGCNA studies, this analysis successfully figured out several mRNAs and lncRNAs correlated to asthma treatment. Functional analysis on these transcripts within the desensitizationtreatment-related module revealed four genes (TRIO, PEA15, KLRK1, NDUFA13) significantly enriched in the positive regulation of apoptotic process. This finding actually was consistent to previous studies for asthma [39][40][41]. Concretely, KLRK1 encodes NKG2D, and aberrant expression of NKG2D ligands has been reported in sites of inflammation and in tissues undergoing autoimmune pathologies, including asthma [44]. Also, Klrk1 −/− mice exhibited a profoundly impaired inflammatory response to allergen challenge compared with klrk1 +/+ mice [36] . The qRT-PCR experiment for KLRK1 based on our samples indicated that the expression of KLRK1 was up-regulated after desensitization treatment, which might alleviate the inflammatory response and respiratory symptoms. This result was consistent to our differential expression analysis, despites the expression difference is not very statistically significant (p value = 0.09676, Figure S7-B). The possible reason may be that we just used three individuals for the qRT-PCR experiments, the sample size might limit the statistical power. The reason why we just used three individuals is that the blood sample was preserved too long despites we kept them refrigerated (− 80°C), the RNAs of some samples are degraded at different degrees. We just found the RNAs of three samples were qualified for qRT-PCR experiments after the quality inspection. Another two transcripts (LINC02145, GUSBP2) for qRT-PCR validation had the same issue. Thereby, it is reasonable that these validations might be more significant if we use ten sample for qRT-PCR experiments. PEA15 has been proved as an important protein that regulates death receptor induced apoptosis and proliferation signaling by binding to FADD and relocating ERK1/2 to the cytosol. It regulates subcellular compartmentalization and activity of phosphoERK, and ERK1 was found to be important for Th2 differentiation and development of experimental asthma [37]. In addition, Sandra Pastorino et al's study showed that ablation of PEA-15 led to accumulation of ERK in the nucleus in PEA-15/ T cells, which would enhance the T-cell proliferation [45]. In asthma, stimuli can promote repair mechanisms to induce cellular proliferation and inhibits apoptosis, and delayed survival of inflammatory cells may in turn aggravate the respiratory symptoms. In the present study, our analysis indicated that the apoptotic process of asthmatic children was upregulated after desensitization treatment, which implies the asthmatic symptoms in those patients might be improved in different degrees. Notably, our analysis figured out a key lncRNA LINC02145 directly targeted on PEA15 (Fig. 7b). Subsequent qRT-PCR on LINC02145 showed the expression of this lncRNA exhibited significant upregulation ( Figure S7-C) after desensitization treatment, which might affect the expression of PEA15 and further regulates death receptor induced apoptosis and proliferation signaling and relieve the respiratory symptoms. Additionally, our analysis indicated that a pseudogene GUSBP2 (also defined as lncRNA in GeneCards Suite [46]) was found to interact with KLRK1, as well as another two immune/inflammatory-related genes, including HMGN1 and WDFY1. HMGN1 encodes the nucleosome-binding protein, which was reported as a potent alarmin that binds TLR4 and induces antigenspecific Th1 immune responses [27]. WDFY1 could mediate Toll-like receptor3/4 (TLR3/4) signaling by recruiting TRIF. TLRs were found to be expressed in the lung tissue and some cells of innate and adaptive immune system. Abnormal expression of TLRs could lead to aberrant expression of many inflammatory and antiinflammatory mediators [36]. These 3 genes thereby have been well demonstrated to be closely associated to immune and inflammatory response in asthmatic status. Regarding GUSBP2, the associations with asthma actually has been reported in a microarray study for asthma, which showed that GUSBP2 was differentially expressed in comparison of asthmatic patients with health controls [47]. However, its detailed mechanisms in asthma was not illuminated. In the present study, our network analysis disclosed GUSBP2 exhibiting correlations to these immune/inflammatory-related genes, implying it might target to these genes to participate the regulation of immune and inflammatory process in asthma. Indeed, GUSBP2 was found as the highest up-regulated fold change of lncRNA in CD (Crohn's disease) patient plasma via microarray screening [48]. CD is known to us as one of inflammatory bowel diseases, which suggests the existence of strong association between GUSBP2 and inflammatory reaction. In the present study, our qRT-PCR experiment on GUSBP2 indicated that the expression of GUSBP2 was up-regulated after desensitization treatment ( Figure S7-A). This finding was also consistent with our differential expression analysis. The result indicated that the expression of GUSBP2 was repressed in asthmatic patients compared to healthy controls (Figure S8). Subsequent network analysis based on PCC method displayed that GUSBP2 might target on KLRK1, whose expression was also already found to be upregulated after desensitization treatment. Notably, the cutoff of correlation value was not set to traditional threshold (0.9~0.95) since the correlation between those transcripts in this study is not as strong as our previous lncRNA studies [38,43]. One of possible reason is that the desensitization treatment might not make a big change of RNAs in peripheral blood of these asthmatics, which further led to weak correlations among those RNAs. This also is the reason that it is not easily to figure out key RNAs related to asthma or treatment simply via differential expression analysis, thereby we applied co-expression network analysis, a  relatively robust and sensitive strategy for our study. Our findings indicated that the desensitization treatment might induce GUSBP2 positively work with KLRK1 to alleviate the inflammatory response and the respiratory symptoms. To validate this further, we analyzed a expression data which derives from GEO database (accession number: GSE2125), which sequenced 45 alveolar macrophages from human subjects (30 non-smoker healthy people and 15 asthmatic patients). The result showed that both genes exhibited a higher expression level in healthy controls ( Figure  S8). Subsequent study on relationship between GUSBP2 and KLRK1 using PCC method indicated that they have a moderate correlation (PCC value = 0.538), which was consistent with It was also implied that GUSBP2 and KLRK1 might play important roles in the occurrence of asthma and might be potential useful RNA molecules for therapy and diagnosis in asthma.

Conclusions
Our network analysis based on the deep RNA sequencing of ten childhood asthmatics before and after desensitization treatment disclosed several lncRNAs exhibiting close correlation to the asthma-related genes, particularly involved in immune, inflammatory response and apoptosis process. This finding provides many promising key noncoding RNA, i.e., LINC01771, LINC02145, GUSBP2, etc. which might be served as novel molecules for observing and supervising desensitization therapy for childhood asthma.

Sampling and deep sequencing
All individuals consented, and the project were approved by the institutional review boards of the Guangzhou Medical University, University of Macau. Ten children diagnosed with asthma in the 1st Affiliated Hospital of Guangzhou Medical University, China were selected for our study. Of the allergic asthma patients had received Dermatophagoides pteronyssinus (Der p) SCIT (before treatment and after half a year). All patients completed questionnaires on their demographics and medical history and fulfilled the Allergic Rhinitis and its Impact on Asthma (ARIA) criteria for allergic rhinitis and/or the Global Initiative for Asthma (GINA) criteria [49] for mild-to-moderate asthma, sIgE to Der p, and Dermatophagoides farinae (Der f). The detailed information for these patients was summarized in Table 4. For each individual, peripheral whole blood was extracted and then the peripheral blood mononuclear cell (PBMC) was separated from it using ficoll-paque. After that, the total

SCIT and drug treatments protocol
The patients were treated with conventional schedule injections of standardized aluminum formulated Der p extract (Alutard, ALK-Abelló, Hørsholm, Demark). The treatment protocol followed the recommended updosing schedule of 17 weeks before reaching a maintenance dose of 100,000 Alutard-SQ ( Figure S6). As-needed use of short-acting bronchodilators (Salbutamol Sulfate Aerosol) for relieving asthma symptoms was also allowed.
For the remaining transcripts which were annotated by aligning against known protein sequences from NCBI nr database, Uniprot database, and the known human mRNAs derived from Ensembl database were excluded. This step aims at eliminating protein-coding sequences as much as possible. The remaining unmapped transcripts were further filtered to remove the sequences with length less than 200 nt and longest ORF longer than 100 residues. Then the second filtering round aimed at removing protein-coding sequences was conducted using Pfamscan [53] and CPC [54]. The remaining transcripts were retained as final dataset of lncRNAs.

Differential expression analysis
All the assembled transcripts for each individual were quantified by featurecounts (v1.5.3) [55]. The expression profile was normalized using median of ratios method by R package DESeq2 [22]. Then the differential expression analysis was performed using pair-wise method of DESeq2 to detect the differentially expressed transcripts (DETs). Functional enrichment analysis for these DETs were conducted based on Metascape serve [56].

Network analysis predicts lncRNA functions
The normalized expression profile was subjected to a R package WGCNA [57] for the weighted coexpression network analysis. Briefly, the function soft-Connectivity from WGCNA was used and the "randomly selected genes" parameter set at 5000, the power parameter pre-calculated by pickSoft-Threshold function of WGCNA. These analyses aim at detection of trait-related modules. Then the selected module with highly correlated RNAs were subjected to Metascape to search for enriched biological processes and pathways. The intramodular connectivity of mRNA and lncRNA in trait-related modules was assessed using Pearson Correlation Coefficient (PCC) method, |cor| > =0.7, p-value<=0.05.

GEO data
The expression profile data for the validation of GUSBP2 and KLRK1 was downloaded from NCBI GEO database, which contains 45 alveolar macrophages from human subjects (30 non-smoker healthy people and 15 asthmatic patients).