- Research article
- Open Access
Transcriptome analysis of peripheral whole blood identifies crucial lncRNAs implicated in childhood asthma
BMC Medical Genomics volume 13, Article number: 136 (2020)
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.
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.
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.
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.
Asthma is a chronic disease of the airways of the lungs, characterized by various symptoms, such as wheezing, coughing, chest tightness, shortness of breath, dyspnea, etc. [1, 2]. It has become one of the most common long-term conditions in children and severely affects the quality of life of children as well as their parents. In Europe, 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 . In China, the prevalence of childhood asthma increased extremely since 1990s, which ranged from 0.93% in 1990 to 1.54% in 2000 . 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 . 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 . 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 . 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 . 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 . 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 . 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 . 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 . The results showed that there was no significant discrepancy in those five AS modes before and after treatment (Fig. 3d).
Differential expression analysis reveals desensitization-treatment-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 . 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 desensitization-treatment-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-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.
Functional analysis of module reveals key genes related to asthma
In order to investigate the potential relationship of transcripts which co-localizes in desensitization-treatment-related 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 (Fig. 5) by four genes: TRIO, PEA15, KLRK1, NDUFA13.
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 . 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 . Furthermore, two genes (G2E3, ZNF19) were found to be differentially methylated between asthmatics and non-asthmatics [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 desensitization-treatment-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 desensitization-treatment-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 MSTRG.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).
With the advancement of high throughput sequencing technology, tens of thousands of lncRNAs have been detected in human. Despite the function of majority of 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 (TH) 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 desensitization-treatment-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 desensitization-treatment-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 . Also, Klrk1−/− mice exhibited a profoundly impaired inflammatory response to allergen challenge compared with klrk1+/+ mice  . 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 . 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 . 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 up-regulation (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 ) 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 antigen-specific Th1 immune responses . 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 anti-inflammatory mediators . 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 . 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 . 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 up-regulated 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.
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  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 RNA was extracted using trizol (invitrogen) method and the RNA quality was assessed using ND-1000 Nanodrop and Agilent 2200 TapeStation. the library for deep RNA sequencing was prepared according to a standard protocol established by RiboBio company in Guangzhou. The rRNAs were removed using Ribo-Zero™ Rrna Removal Reagent (Human/Mouse/Rat)-Illumina and library was prepared using NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina, USA NEB company. Then the paired-end 150 bp RNAs sequencing was performed using Illumina HiseqX-ten.
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.
Basic bioinformatic analysis for identification of lncRNAs
The raw RNA-seq data were filtered using Trimmomatic v0.36  (ILLUMINACLIP: TruSeq3-PE.fa:2:30:10:8:true SLIDINGWINDOW:4:15 LEADING:3 TRAILING:3 MINLEN:50) to remove low quality reads and adaptors. Then the clean reads were aligned to Ensembl hg38 human genome via STAR (v020201) . Then with the help of program StringTie (v1.3.3b) , the output bam files were utilized to re-construct the new transcriptome for this study.
The assembled transcripts were subjected to stringent stepwise filtering pipeline to identify the high-confidence dataset of lncRNAs, which has been applied in our previous studies. Briefly, the known lncRNAs were initially extracted based on the “biotype_transcript” of reference gtf file of Homo sapiens from all the assembled transcripts via an in-house perl script. “biotype_transcript” tagged as “Long_non-coding_RNA”, “Non_coding”, “3prime_overlapping_ncRNA”, “Antisense”, “lincRNA”, “Retained_intron”, “Sense_intronic”, “Sense_overlapping”, “macro_lncRNA”, “bidirectional_promoter_lncRNA” was detected as lncRNAs.
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  and CPC . 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) . The expression profile was normalized using median of ratios method by R package DESeq2 . 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 .
Network analysis predicts lncRNA functions
The normalized expression profile was subjected to a R package WGCNA  for the weighted co-expression network analysis. Briefly, the function softConnectivity 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.
Validation by qRT-PCR
To validate the RNA expression level for the lncRNAs or mRNAs, we determined the expression levels of several selected lncRNA or mRNA transcripts by RT-qPCR in two different statuses (before and after treatment), including KLRK1, LINC02145 and GUSBP2. Primers were designed and synthesized by RiboBio company in Guangzhou. Upstream primer of KLRK1 is: 5′- CCGACACAAAGTCCCACACTC − 3′, downstream primer: 5′- AGATGCTTGCCTAAACGCCT − 3′. Upstream primer of LINC02145 is: 5′-AAGTACTGTTCTGCCTTCCACAC-3′, downstream primer: 5′-TGCCTTGGAAAGGTAGCGAG-3′. Upstream primer of GUSBP2 is: 5′- GAGCAGTACCATCTGGGTCTGG-3′, downstream primer: 5′- ACTTCATCTTGGATTTCCAGCCT-3′. Prepare 10 μl SYBR Green qRT-PCR kit (Soochow GenePharma Co. Ltd.) reaction system: cDNA (500 ng/20 μl): 1 μl, forward primer (10 pmol/L): 0.4 μl, reverse primer (10 pmol/L): 0.4 μl, 2 × SYBR Green PCR Master Mix 5 μl, RNase-free H2O 3.2 μl. Reaction condition: 95 °C 10 min heating, 95 °C 15 s denaturing, 60 °C 30 s annealing, 70 °C 30 s extension, 40 cycles. All reactions had 3 duplicates, CCAT1, E-cadherin and N-cadherin expression in colorectal tissue was calculated with 2-ΔΔCt method. One-side T-test was used to test whether there is a difference in the expression level between groups before and after desensitization treatment. These analyses were achieved on the open source statistical programing language R. Concretely, the function shapiro.test() was used to test whether two groups before and after desensitization treatment are normal distribution. The function var.test() was used to test whether two groups are equal variance. The function t.test() was used to perform paired t-test for two groups.
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).
Availability of data and materials
All the raw data used in this study is deposited at DDBJ/EMBL/GenBank under the SRA accession PRJNA552157, associated with the BioSample SAMN12179172-SAMN12179191.
Long noncoding RNA
Airway smooth muscle
T helper cell 1
Open reading frame
- TH cell:
Helper T cell
Methyl CpG binding protein 1
Pearson Correlation Coefficient
Peripheral blood mononuclear cell
Differentially expressed transcripts
Gene Ontology; nr database
Buelo A, McLean S, Julious S, Flores-Kim J, Bush A, Henderson J, Paton JY, Sheikh A, Shields M, Pinnock H, et al. At-risk children with asthma (ARC): a systematic review. Thorax. 2018;73(9):813–24.
van Aalderen WM. Childhood asthma: diagnosis and treatment. Scientifica. 2012;2012:674204.
Lin R, Guan R, Liu X, Zhao B, Guan J, Lu L. Significant rise of the prevalence and clinical features of childhood asthma in Qingdao China: cluster sampling investigation of 10,082 children. BMC Public Health. 2014;14:1002.
Konradsen JR, Caffrey Osvald E, Hedlin G. Update on the current methods for the diagnosis and treatment of severe childhood asthma. Expert Rev Respir Med. 2015;9(6):769–77.
Singh V. What is new in the management of childhood asthma? Indian J Pediatr. 2008;75(8):845–53.
Nordlund B, Melen E, Schultz ES, Gronlund H, Hedlin G, Kull I. Prevalence of severe childhood asthma according to the WHO. Respir Med. 2014;108(8):1234–7.
Arrieta MC, Stiemsma LT, Dimitriu PA, Thorson L, Russell S, Yurist-Doutsch S, Kuzeljevic B, Gold MJ, Britton HM, Lefebvre DL, et al. Early infancy microbial and metabolic alterations affect risk of childhood asthma. Sci Transl Med. 2015;7(307):307ra152.
Chen Z, Salam MT, Alderete TL, Habre R, Bastain TM, Berhane K, Gilliland FD. Effects of childhood asthma on the development of obesity among school-aged children. Am J Respir Crit Care Med. 2017;195(9):1181–8.
McConnell R, Berhane K, Yao L, Jerrett M, Lurmann F, Gilliland F, Künzli N, Gauderman J, Avol E, Thomas D. Traffic, susceptibility, and childhood asthma. Environ Health Perspect. 2006;114(5):766–72.
Weiss ST, Silverman EK. Pro: genome-wide association studies (GWAS) in asthma. Am J Respir Crit Care Med. 2011;184(6):631–3.
Feil R, Fraga MF. Epigenetics and the environment: emerging patterns and implications. Nat Rev Genet. 2012;13(2):97–109.
Aguilera O, Fernandez AF, Munoz A, Fraga MF. Epigenetics and environment: a complex relationship. J Appl Physiol. 2010;109(1):243–51.
Austin PJ, Tsitsiou E, Boardman C, Jones SW, Lindsay MA, Adcock IM, Chung KF, Perry MM. Transcriptional profiling identifies the long noncoding RNA plasmacytoma variant translocation (PVT1) as a novel regulator of the asthmatic phenotype in human airway smooth muscle. J Allergy Clin Immunol. 2017;139(3):780–9.
Khalil AM, Guttman M, Huarte M, Garber M, Raj A, Rivea Morales D, Thomas K, Presser A, Bernstein BE, van Oudenaarden A, et al. Many human large intergenic noncoding RNAs associate with chromatin-modifying complexes and affect gene expression. Proc Natl Acad Sci U S A. 2009;106(28):11667–72.
Rinn JL, Chang HY. Genome regulation by long noncoding RNAs. Annu Rev Biochem. 2012;81:145–66.
Gomez JA, Wapinski OL, Yang YW, Bureau JF, Gopinath S, Monack DM, Chang HY, Brahic M, Kirkegaard K. The NeST long ncRNA controls microbial susceptibility and epigenetic activation of the interferon-gamma locus. Cell. 2013;152(4):743–54.
Gupta RA, Shah N, Wang KC, Kim J, Horlings HM, Wong DJ, Tsai MC, Hung T, Argani P, Rinn JL, et al. Long non-coding RNA HOTAIR reprograms chromatin state to promote cancer metastasis. Nature. 2010;464(7291):1071–6.
Hu G, Tang Q, Sharma S, Yu F, Escobar TM, Muljo SA, Zhu J, Zhao K. Expression and regulation of intergenic long noncoding RNAs during T cell development and differentiation. Nat Immunol. 2013;14(11):1190–8.
Collier SP, Collins PL, Williams CL, Boothby MR, Aune TM. Cutting edge: influence of Tmevpg1, a long intergenic noncoding RNA, on the expression of Ifng by Th1 cells. J Immunol. 2012;189(5):2084–8.
Zhang H, Nestor CE, Zhao S, Lentini A, Bohle B, Benson M, Wang H. Profiling of human CD4+ T-cell subsets identifies the TH2-specific noncoding RNA GATA3-AS1. J Allergy Clin Immunol. 2013;132(4):1005–8.
Foissac S, Sammeth M. ASTALAVISTA: dynamic and flexible analysis of alternative splicing events in custom gene datasets. Nucleic Acids Res. 2007;35(Web Server issue):W297–9.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.
Waterfield M, Khan IS, Cortez JT, Fan U, Metzger T, Greer A, Fasano K, Martinez-Llordella M, Pollack JL, Erle DJ. The transcriptional regulator Aire coopts the repressive ATF7ip-MBD1 complex for the induction of immunotolerance. Nat Immunol. 2014;15(3):258.
Zhang Y. Potential therapeutic targets from genetic and epigenetic approaches for asthma; 2016.
Liang L, Willis-Owen SAG, Laprise C, Wong KCC, Davies GA, Hudson TJ, Binia A, Hopkin JM, Yang IV, Grundberg E, et al. An epigenome-wide association study of total serum immunoglobulin E concentration. Nature. 2015;520(7549):670–4.
Poppinga WJ, Munoz-Llancao P, Gonzalez-Billault C, Schmidt M. A-kinase anchoring proteins: cAMP compartmentalization in neurodegenerative and obstructive pulmonary diseases. Br J Pharmacol. 2014;171(24):5603–23.
Wei F, Yang D, Tewary P, Li Y, Li S, Chen X, Howard OM, Bustin M, Oppenheim JJ. The Alarmin HMGN1 contributes to antitumor immunity and is a potent immunoadjuvant. Cancer Res. 2014;74(21):5989–98.
Papaioannou AI, Spathis A, Kostikas K, Karakitsos P, Papiris S, Rossios C. The role of endosomal toll-like receptors in asthma. Eur J Pharmacol. 2017;808:14–20.
Zhang Q, Fu XL, Qian FH, Cao Q, Mao ZD, Bai JL, Du Q, Shi Y. Polymorphisms in toll-like receptor 3 are associated with asthma-related phenotypes in the Chinese Han patients. Int J Immunogenet. 2016;43(6):383–90.
Wagener AH. Biomarker discovery for asthma phenotyping: from gene expression to the clinic: 9789461697783; 2016.
Goplen N, Karim Z, Guo L, Zhuang Y, Huang H, Gorska MM, Gelfand E, Pagés G, Pouysségur J, Alam R. ERK1 is important for Th2 differentiation and development of experimental asthma. FASEB J. 2012;26(5):1934–45.
Jevnikar Z, Ostling J, Ax E, Calven J, Thorn K, Israelsson E, Oberg L, Singhania A, Lau LCK, Wilson SJ, et al. Epithelial IL-6 trans-signaling defines a new asthma phenotype with increased airway inflammation. J Allergy Clin Immunol. 2019;143(2):577–90.
Ho E, Dagnino L. Emerging role of ILK and ELMO2 in the integration of adhesion and migration pathways. Cell Adhes Migr. 2012;6(3):168–72.
Murphy TM, Wong CC, Arseneault L, Burrage J, Macdonald R, Hannon E, Fisher HL, Ambler A, Moffitt TE, Caspi A, et al. Methylomic markers of persistent childhood asthma: a longitudinal study of asthma-discordant monozygotic twins. Clin Epigenetics. 2015;7:130.
Rastogi D, Suzuki M, Greally JM. Differential epigenome-wide DNA methylation patterns in childhood obesity-associated asthma. Sci Rep. 2013;3:2164.
Farhadi N, Lambert L, Triulzi C, Openshaw PJ, Guerra N, Culley FJ. Natural killer cell NKG2D and granzyme B are critical for allergic pulmonary inflammation. J Allergy Clin Immunol. 2014;133(3):827–35 e823.
Kerbrat S, Vingert B, Junier MP, Castellano F, Renault-Mihara F, Dos Reis TS, Surenaud M, Noizat-Pirenne F, Boczkowski J, Guellaen G, et al. Absence of the adaptor protein PEA-15 is associated with altered pattern of Th cytokines production by activated CD4+ T lymphocytes in vitro, and defective red blood cell Alloimmune response in vivo. PLoS One. 2015;10(8):e0136885.
Huang C, Leng D, Sun S, Zhang XD. Re-analysis of the coral Acropora digitifera transcriptome reveals a complex lncRNAs-mRNAs interaction network implicated in Symbiodinium infection. BMC Genomics. 2019;20(1):48.
Todo-Bom A, Mota Pinto A, Alves V, Vale Pereira S, Santos Rosa M. Apoptosis and asthma in the elderly. J Investig Allergol Clin Immunol. 2007;17(2):107–12.
Zhou C, Yin G, Liu J, Liu X, Zhao S. Epithelial apoptosis and loss in airways of children with asthma. J Asthma. 2011;48(4):358–65.
O’Sullivan MP, Tyner JW, Holtzman MJ. Apoptosis in the airways: another balancing act in the epithelial program. Am J Respir Cell Mol Biol. 2003;29(1):3–7.
Huang C, Morlighem JRL, Cai J, Liao Q, Perez CD, Gomes PB, Guo M, Radis-Baptista G, Lee SM. Identification of long non-coding RNAs in two anthozoan species and their possible implications for coral bleaching. Sci Rep. 2017;7(1):5333.
Huang C, Leng D, Lei KC, Sun S, Zhang XD. Transcriptome analysis reveals lncRNA-mediated complex regulatory network response to DNA damage in the liver tissue of Rattus norvegicus. J Cell Physiol. 2019;234(12):23216–31.
Lanier LL. NKG2D receptor and its ligands in host defense. Cancer Immunol Res. 2015;3(6):575–82.
Pastorino S, Renganathan H, Caliva MJ, Filbert EL, Opoku-Ansah J, Sulzmaier FJ, Gawecka JE, Werlen G, Shaw AS, Ramos JW. The death effector domain protein PEA-15 negatively regulates T-cell receptor signaling. FASEB J. 2010;24(8):2818–28.
Stelzer G, Rosen N, Plaschkes I, Zimmerman S, Twik M, Fishilevich S, Stein TI, Nudel R, Lieder I, Mazor Y. The GeneCards suite: from gene data mining to disease genome sequence analyses. Curr Protoc Bioinformatics. 2016;54(1):1.30. 31–31.30. 33.
Orsmark-Pietras C, James A, Konradsen JR, Nordlund B, Söderhäll C, Pulkkinen V, Pedroletti C, Daham K, Kupczyk M, Dahlén B. Transcriptome analysis reveals upregulation of bitter taste receptors in severe asthmatics. Eur Respir J. 2013;42(1):65–78.
Chen D, Liu J, Zhao H-Y, Chen Y-P, Xiang Z, Jin X. Plasma long noncoding RNA expression profile identified by microarray in patients with Crohn’s disease. World J Gastroenterol. 2016;22(19):4716.
Program GIfA. Global Strategy for Asthma Management and Prevention: 2012 Update. In: Global Initiative for Asthma Vancouver, USA; 2017.
Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30(15):2114–20.
Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, Batut P, Chaisson M, Gingeras TR. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29(1):15–21.
Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. 2015;33(3):290–5.
Chojnacki S, Cowley A, Lee J, Foix A, Lopez R. Programmatic access to bioinformatics tools from EMBL-EBI update: 2017. Nucleic Acids Res. 2017;45(W1):W550–3.
Kong L, Zhang Y, Ye Z-Q, Liu X-Q, Zhao S-Q, Wei L, Gao G. CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Res. 2007;35(suppl_2):W345–9.
Liao Y, Smyth GK, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 2014;30(7):923–30.
Zhou Y, Zhou B, Pache L, Chang M, Khodabakhshi AH, Tanaseichuk O, Benner C, Chanda SK. Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat Commun. 2019;10(1):1523.
Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559.
Thanks to Shixue Sun and Chang Chen for the suggestions of manuscript revision.
This work was funded by University of Macau (grant numbers: FHS-CRDA-029-002-2017, EF005/FHS-ZXH/2018/GSTIC, MYRG2018–00071-FHS and SRG2016–00083-FHS), National Natural Science Foundation of China (Project No.: 81871736, 81601394, 81572063), Bureau of traditional Chinese Medicine Scientific Research Project of Guangdong (Project No.: 20192048), Research Project of First Affiliated Hospital of Guangzhou Medical University (Project No.: ZH201915) and by The Science and Technology Development Fund, Macau SAR (File no. 0004/2019/AFJ). Individuals working for the funding body played a role in the data analysis and in manuscript writing, as outlined in the Author’s contributions.
Ethics approval and consent to participate
This study was approved by the Medical Ethics Committee of First Affiliated Hospital of Guangzhou Medical University (ethics approval no. gyfyy-2016-73). All experiments were performed in accordance with relevant guidelines and regulations of the Ethics Committee of First Affiliated Hospital of Guangzhou Medical University. The informed consent for participation in the study was obtained from their parents.
Consent for publication
All the parents and legal guardians consented to the publication of their children of clinical and sequencing data.
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Additional file 1: Table S1.
Basic statistics of deep RNA sequencing data before and after processing.
Additional file 2: Table S2.
Basic statistics of clean reads alignment against to the Homo sapiens genome sequence.
Additional file 3: Table S3.
Basic statistics of assembly results of transcriptome in Homo sapiens.
Additional file 4: Table S4.
Location of novel identified lncRNAs.
Additional file 5: Figure S1.
Sample clustering to detect outliers. All the samples were in the clusters, all samples have passed the cuts.
Additional file 6: Figure S2.
Analysis of network topology for various soft-thresholding powers. The left panel indicates the scale-free fit index (y-axis) as a function of the soft-thresholding power (x-axis). The right panel shows the mean connectivity (degree, y-axis) as a function of the soft-thresholding power (x-axis).
Additional file 7: Figure S3.
Clustering dendrograms of transcripts, with dissimilarity based on topological overlap, together with assigned module colors.
Additional file 8: Figure S4.
Visualizing the gene network using a heatmap plot. Light color represents low overlap and progressively darker red color represents higher overlap.
Additional file 9: Figure S5.
Scatterplots of Gene Significance (GS) for recurrence vs Module Membership (MM) in the desensitization treatment-related module.
Additional file 10: Figure S6.
The treatment protocol applied in the present study.
Additional file 11: Figure S7.
QRT-PCR validation for three selected key lnRNAs or mRNAs. P value was calculated by paired T-test.
Additional file 12: Figure S8.
Expression comparison between GUSBP2 and KLRK1 in GEO GSE2125. P value was calculated by unpaired T-test.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Zheng, P., Huang, C., Leng, D. et al. Transcriptome analysis of peripheral whole blood identifies crucial lncRNAs implicated in childhood asthma. BMC Med Genomics 13, 136 (2020). https://doi.org/10.1186/s12920-020-00785-y
- Childhood asthma
- Deep RNA-sequencing
- Long non-coding RNAs