The mRNA and miRNA expression profiles of asthma patients were obtained from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) utilizing the getGEO function of the GEOquery package (version 2.58.0) in R software (version 4.0.3) (Barrett et al., 2012). The following terms were used to search the microarray studies: “Asthma”, AND “epithelial cells”, OR “bronchial brushings”, AND “mRNA expression”, OR “miRNA expression”, AND “Homo sapiens”, AND “Microarray”. The following eligibility criteria were used to include or exclude datasets and samples: 1) A minimum of 10 subjects in the dataset, containing at least four asthma patients and four controls; 2) Bronchial brushings were used for microarray analysis; 3) All subjects were older than 18 years old and were steroid-naïve; 4) Raw data were available in GEO.

The expression matrix file and related platform annotation document of each microarray dataset were downloaded from the GEO database to convert the names of microarray probes to the gene symbols. Probes matching with multiple gene symbols were eliminated and the mean values were calculated for gene symbols corresponding to multiple probes. After excluding samples with strong heterogeneity according to the hclust method, the differentially expressed genes (DEGs) and the differentially expressed miRNAs (DEMis) between asthma and control samples in each dataset were identified by using the “limma” (linear models for microarray data) package (version 3.46.0) (Ritchie et al., 2015) in R software with the cut-off criteria of |log2 fold change| > 0.5 and p-value < 0.05.

To integrate the outcomes of multiple microarray datasets and minimize the bias and inconsistencies, the RRA method was adopted to identify the robust DEGs. The up- and down-regulated DEGs were firstly ranked by expression fold changes in each dataset and were subsequently analyzed using the “Robust Rank Aggregation” R package (version 1.1) (Kolde et al., 2012). The score in the RRA result indicated the ranking degree of each gene in the final gene list. Genes with score < 0.05 and |log2 fold change| > 0.5 were considered as the significant robust DEGs.

The “RCircos” R package (1.2.1) was utilized to visualize the expression patterns and the chromosomal positions of all robust DEGs.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were used to investigate the biological process (BP), the cellular component (CC), the molecular function (MF), and the involved pathways of selected molecules, which were performed with the “clusterprofiler” R package (version 3.18.1) (Yu et al., 2012). The GO terms and KEGG pathways with p-value < 0.05 were considered statistically significant and further visualized via the “ggplot2” R package (version 3.3.3). The KEGG enrichment analysis of differentially expressed miRNAs was conducted using the miRPathDB v2.0 (https://mpd.bioinf.uni-sb.de/) database (Kehl et al., 2020).

DisGeNET (http://www.disgenet.org) database is one of the largest available platforms of human disease-associated genes and variants through manually integrating the scientific literature (Piñero et al., 2017). For a given gene list, DisGeNET database can identify significantly correlated diseases.

All previously identified robust DEGs were uploaded to the STRING (version 11.0) database (https://www.string-db.org/) to construct the protein-protein interaction (PPI) network (Szklarczyk et al., 2021). Confidence > 0.4 was set as the screening criteria. The PPI network was subsequently reconstructed and visualized via the Cytoscape (version 3.8.2) (http://cytoscape.org/) software (Su et al., 2014). In the Cytoscape plot, each node represented a gene/protein/miRNA/circRNA, while the edge between nodes represented the interactions of molecules. The molecular complex detection (MCODE) plugin of the Cytoscape software was used to screen out significant clusters in the PPI network.

CytoHubba is another plugin of the Cytoscape software for ranking nodes in a network, which provides eleven topological analysis methods and six centralities to identify hub genes based on shortest paths, including Maximal Clique Centrality (MCC), Density of Maximum Neighborhood Component (DMNC), Maximum Neighborhood Component (MNC), Degree, Edge Percolated Component (EPC), BottleNeck, EcCentricity, Closeness, Radiality, and Betweenness (Chin et al., 2014). Considering the potential heterogeneity of the biological network, the RRA method was adopted to integrate the results of different analysis methods for catching essential proteins.

GeneMANIA (http://www.genemania.org) database was used to construct the PPI network and explore the putative functions of up-loaded genes (Warde-Farley et al., 2010). For a given query list, GeneMANIA analyzes target genes with functionally similar genes together to obtain regulatory networks.

The Targetscan (http://www.targetscan.org/vert_72/) (Garcia et al., 2011), miRDB (http://mirdb.org/) (Chen and Wang, 2020), and miRWalk (http://mirwalk.umm.uni-heidelberg.de/) (Dweep et al., 2011) databases were used to predict the corresponding miRNAs of all 19 hub genes. The overlapping results of three databases were then intersected with the differentially expressed miRNAs of GSE142237. The ENCORI (http://starbase.sysu.edu.cn/index.php) database (Li et al., 2014) was used to predict the upstream circRNAs of the selected miRNA-mRNA pairs. The final ceRNA network was further processed with the Cytoscape software (version 3.8.2).

The differential analysis was conducted by the “limma” package (version 3.46.0) in R version 4.0.3. Heatmap was used to reveal the logarithmic fold changes of robust DEGs in the RRA analysis. p < 0.05 was considered statistically significant.

Five mRNA microarray datasets (GSE4302, GSE43696, GSE63142, GSE67472, and GSE41861) and one miRNA microarray dataset (GSE142237) derived from bronchial epithelial brushings were obtained from the GEO database. There were a total of 272 steroid-naïve asthma patients and 165 healthy controls in the five mRNA microarray datasets. The miRNA microarray dataset (GSE142237) included a total of eight asthma patients and four healthy controls. Only asthma patients without any steroid treatments were included for further analysis.

The workflow of the study was shown in Figure 1. Detailed information on the datasets mentioned above was shown in Table 1.

After base two logarithm conversions, all five mRNA microarray datasets were further analyzed using the “limma” package of the R software to identify DEGs in primary airway epithelial cells between steroid-naïve asthma patients and healthy controls. The cut-off criteria were as follows: p-value < 0.05 and |log2 fold change| > 0.5. The distribution of DEGs in five datasets was illustrated by the volcano plots in Figure 2. There were 851 DEGs (737 upregulated and 114 downregulated), 1218 DEGs (606 upregulated and 612 downregulated), 783 DEGs (345 upregulated and 438 downregulated), 1084 DEGs (650 upregulated and 434 downregulated), and 875 DEGs (348 upregulated and 527 downregulated) in GSE4302, GSE43696, GSE63142, GSE67472, and GSE41861, respectively.

The RRA analysis was conducted to integrate five mRNA microarray datasets with minimal bias and a total of 127 integrated DEGs were identified, including 83 upregulated and 44 downregulated robust DEGs. Lower RRA scores in the results indicated higher gene ranks and the top 20 genes were shown in the lollipop chart (Figure 3A), only five of which were downregulated. The top 10 ranked genes in asthma were composed of eight upregulated genes [including CLCA1 (score = 4.99E-15), CPA3 (score = 1.21E-12), CST1 (score = 1.56E-11), PRR4 (score = 2.51E-11), SERPINB2 (score = 1.08E-08), POSTN (score = 3.57E-08), ITLN1 (score = 7.85E-08), UPK1B (score = 1.11E-07)] and two downregulated genes [including BPIFA1 (score = 3.52E-12), FHOD3 (score = 3.90E-08)]. Of note, CLCA1, SERPINB2, and POSTN, as the IL-13-responsive genes and airway epithelial gene signatures for type 2 status, were all at a higher rank in five microarray datasets. This suggested the significant roles of the three genes in asthma pathogenesis and the reliability of the RRA results (Woodruff et al., 2007; Woodruff et al., 2009). The top 10 up- and down-regulated DEGs were drawn on the heatmap (Figure 3B). The expression patterns and the chromosomal positions of all robust DEGs across the five microarray datasets were shown in Figure 3C. These DEGs were distributed in all chromosomes except chromosome Y. Chromosome three contained the most DEGs.

GO annotation and KEGG pathway enrichment analyses of robust DEGs were performed to explore the biological classifications in asthma utilizing the clusterProfiler R package (Figures 4A–C). Three categories including biological process, cellular component, and molecular function were analyzed within up- and down-regulated DEGs, respectively. For 83 upregulated genes, the outcomes showed that changes in the biological process were mainly enriched in regulation of endopeptidase and peptidase activities, followed by response to tissue homeostasis, negative regulation of proteolysis. The significantly enriched entries for the cellular component part were collagen-containing extracellular matrix, cornified envelope, and external side of plasma membrane. Furthermore, endopeptidase and peptidase regulator activities, enzyme inhibitor activity, and cysteine-type endopeptidase inhibitor activity accounted for the majority of the molecular function terms (Figure 4A). In terms of 44 downregulated genes, the significantly enriched biological process terms were humoral immune response, response to drug, and pattern specification process. In the cellular component part, the downregulated genes were particularly enriched in tight junction, brush border membrane, and Z disc. Meanwhile, endopeptidase and peptidase regulator activities, enzyme inhibitor activity, and heme binding were mainly enriched in the molecular function group (Figure 4B). In addition, integrated DEGs were mainly involved in salivary secretion, metabolism of xenobiotics by cytochrome P450, IL-17 signaling pathway, and leukocyte transendothelial migration in KEGG pathway analysis (Figure 4C). The DisGeNET database was further used to identify DEGs related diseases. As shown in Figure 4D, the result indicated that robust DEGs participated in the progression of various diseases, such as Nasal Polyps, Allergic rhinitis disorder, Allergic asthma, and Atopic Dermatitis, which were all related to allergic reactions and chronic inflammation (Figure 4D). Taken together, the above results indicated that the robust DEGs were mostly associated with asthma-related functions.

In order to explore the potential protein-protein interactions in asthma, all 127 robust DEGs were uploaded to the STRING database for further analysis (http://string.embl.de/). After hiding the disconnected nodes, the Cytoscape software was adopted to visualize the network (Figure 5A). As shown in the final network, 77 nodes and 114 edges were obtained, including 57 upregulated and 20 downregulated genes. Three key clusters were identified from the whole network using the MCODE plugin (Figures 5B–D). GO enrichment analyses showed that the significantly enriched biological process terms of three clusters were regulation of myeloid leukocyte mediated immunity, T cell activation, and antibacterial humoral response, respectively (Figure 5E).

Hub genes were subsequently screened out utilizing the cytoHubba plugin, which investigates the most important nodes in the PPI network with several topological analysis algorithms. In order to improve the positive rate of hub gene identification, the RRA method was used to integrate the top 50 ranked genes calculated by ten different topological analysis algorithms (MCC, DMNC, MNC, Degree, EPC, BottleNeck, EcCentricity, Closeness, Radiality, and Betweenness) and a total of 19 genes were obtained accordingly (Figure 5F). The detailed information of all 19 genes was listed in Table 2, including full names, scores of the RRA method, direction, and primary functions. The upset diagram of the top 50 ranked genes from the ten algorithms was shown in Supplementary Figure 1. The GeneMANIA database was used to construct the regulatory network among these 19 hub genes with functionally similar genes and the result showed that these genes may be involved in the following functions: tissue homeostasis, secretory granule, serine-type peptidase and endopeptidase activities, serine hydrolase activity, regulation of protein processing, and humoral immune response (Figure 6).

All 19 hub genes, including 14 upregulated genes (CST4, CTSG, CLCA1, CSTA, CPA3, KIT, SERPINB2, GGH, MUC5AC, POSTN, ADRA2A, TPSAB1, CD69, and AZGP1) and five downregulated genes (LTF, C3, MUC5B, BPIFA1, and SCGB1A1), were used for the circRNA-miRNA-mRNA network construction. MiRNA-mRNA analyses were performed with Targetscan, miRDB, and miRWalk databases. The intersection of miRNA results predicted by three databases was selected as the prediction result. Furthermore, GSE142237 was adopted to validate the prediction result. The volcano plot showed the distribution of DEMis of GSE142237 and there were 522 DEMis (184 upregulated and 338 downregulated) in asthma (Figure 7A). Generally, miRNA has a negative regulatory relationship with its targeted mRNA (Yekta et al., 2004; Lewis et al., 2005; Gebert and MacRae, 2019; Goodall and Wickramasinghe, 2021). Therefore, the miRNAs targeting upregulated hub genes further intersected with downregulated DEMis of GSE142237, while the miRNAs targeting downregulated hub genes intersected with upregulated DEMis. As shown in the Venn diagrams in Figures 7B,C and Table 3, there were 45 miRNA-mRNA relationships in total. The KEGG analyses of seven upregulated and 34 downregulated miRNAs were further conducted via the miRPathDB database (Supplementary Figure 2).

Previous reports have demonstrated that circRNA could function as miRNA sponge to prevent co-expressed mRNA from miRNA-mediated degradation (Piwecka et al., 2017; Yu et al., 2017; Kleaveland et al., 2018; Kristensen et al., 2018; Wang et al., 2019a; Wong et al., 2020). Therefore, the target mRNA has the same expression pattern as circRNA. Consequently, 45 miRNA-mRNA relationships were used to construct the potential circRNA-miRNA-mRNA network. The corresponding circRNAs were predicted with the ENCORI database, while hsa-miR-1293-CST4 and hsa-miR-3664-5p-CTSG failed to obtain targeted circRNAs. After cross-linking, only circRNAs that were able to regulate all remaining mRNAs by targeting corresponding miRNAs were retained to build the network. As was shown in Figure 7D, a ceRNA network, with three circRNAs, 27 miRNAs, and 12 mRNAs, was finally constructed. The detailed information of three circRNAs was listed in Supplementary Table 1. The top 10 ranked ceRNAs identified by MCC algorithm were shown in Figure 7E, including hsa_circ_0001585, hsa_circ_0078031, hsa_circ_0000552, hsa-miR-30a-3p, hsa-miR-30d-3p, KIT, CD69, ADRA2A, BPIFA1, and GGH.