We downloaded the asthma microarray data set GSE89809, including airway epithelial samples from 27 patients with mild to moderate asthma and 11 patients with severe asthma (14), from the GEO (https://www.ncbi.nlm.nih.gov/geo/) database from GPL13158 (Affymetrix HT HG-U133+ PM Array Plate). Raw data were background corrected and normalized by the RMA algorithm. Differentially expressed genes (DEGs) were screened using the “limma” software package (15), and P < 0.05 and |log2FC| > 0.585 were considered to indicate significant differences. ARGs were extracted from the Human Autophagy Database (http://www.autophagy.lu/index.html) and the GO_AUTOPHAGY gene set in the Gene Set Enrichment Analysis (GSEA) website (http://software.broadinstitute.org/gsea/index.jsp). In accordance with Ref., these two gene sets were merged into one ARG set (16).

The unsupervised cluster analysis was performed to identify different subtypes on the basis of the differences in ARG expression in asthma (17, 18). A consensus clustering algorithm was used to evaluate cluster numbers and robustness. The R package “ConsensuClusterPlus” implemented the above steps for 1000 iterations to guarantee the robustness of classification (19). The protein–protein interaction (PPI) network was obtained from the STRING database (https://string-db.org/), and the MCODE was performed using the Metascape (http://metascape.org/gp/index.html#/main/step1) platform. We used the “clusterProfiler” and “ggplot2” packages to perform Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses, respectively, on DEGs (20, 21). GSEA was performed on the gene expression matrix through the “clusterProfiler” package, and the “c2.cp.kegg.v7.0.symbols.gmt” was selected as the reference gene set. Genes associated with specific subtypes were identified by the weighted gene co-expression network analysis (WGCNA) by using the R package “WGCNA” (22). The topological overlap matrix (TOM) and the corresponding dissimilarity (1 − TOM) were transformed from the adjacency matrix. A hierarchical clustering dendrogram was further built, and similar gene expressions were divided into different modules. Finally, the expression profiles of each module were summarized by the module eigengene (ME), and the correlation between the ME and clinical features was determined.

We used the least absolute shrinkage and selection operator (LASSO) logistic regression and the support vector machine–recursive feature elimination (SVM-RFE) to perform feature selection and screen diagnostic markers for asthma (23, 24). The LASSO algorithm was applied with the “glmnet” package (25). The SVM module was established to further identify the diagnostic value of these biomarkers in asthma by the “e1071” package (26). Subsequently, the degree of immune cell infiltration in the airway epithelium was assessed. The CIBERSORT algorithm is an excellent tool used to calculate the abundance of specific cells in the mixture matrix (27). The “Corrplot” software package was further used to draw relevant heatmaps and show the correlation of 22 types of infiltrating immune cells (28), and the differences in infiltrating immune cells between groups were illustrated by the “ggplot2” package violin plot. Finally, the Spearman correlation analysis of diagnostic markers and infiltrating immune cells was performed using the “ggstatsplot” package, and results were visualized using the “ggplot2” package (29).

Univariate Cox regression models were used to select genes associated with asthma control. P < 0.05 was considered statistically significant. The LASSO regression was used to screen out the optimal gene combination and construct the risk characteristics, and the regression coefficient was linearly combined with the gene expression level to establish the risk characteristics. The risk score was calculated as follows: Risk score = (exprgene1 × Coefgene1) + (exprgene2 × Coefgene2) + … + (exprgenen × Coefgenen) (30). Patients with asthma were divided into low- and high-risk groups in accordance with the median risk score. The Kruskal test was used to compare the infiltrating immune cell abundance score and immune response score in the two risk models.

The transcriptional regulatory networks of prognosis-related genes were predicted from the ChEA3 database (https://maayanlab.cloud/chea3/). The database integrated ENCODE; ReMap; some independently published CHIP-seq data; and transcription factor co-expression data from GTEx, TCGA, and ARCHS4 RNA-seq data (31). The target miRNA was first predicted using the starBase (http://starbase.sysu.edu.cn/) database to predict the regulation of target genes by noncoding RNAs, and prediction results included the analysis of RNA22, miRanda, and TargetScan (32). Then, the target lncRNA of miRNA was predicted using the miRNet2.0 database (www.mirnet.ca/miRNet/home.xhtml) and starBase. Finally, the ceRNA network was established (33, 34).

The Broad Institutes Connectivity Map (cMAP) database (https://portals.broadinstitute.org/cmap) was used for the identification of small candidate molecules related to asthma (35). For the identification of small candidate chemical molecules, DEGs (log2FC > 0.585) were introduced into the cMAP database for GSEA. The PubChem (https://pubchem.ncbi.nml.gov) was used for the extraction of detailed information and 3D confirmation of the established small molecules. Then, molecular docking experiments were carried out, and the crystal structures of key targets were retrieved from the human (Human) database in the PDB (https://www.rcsb.org/) database. AutoDockTool software was used to dewater and hydrogenate the receptor protein, and Discovery Studio 4.5.0 software was used to perform molecular docking between small drug molecules and target proteins. A negative value indicates free binding. In this study, binding energy value less than − 5 kcal/mol (20.9 kj/mol) was set as significant binding.

BALB/c mice (female, 17–20 g, 6–8 weeks old) were obtained from Beijing Vital River Laboratory Animal Technology Co., Ltd in China. All animal studies were conducted in accordance with the institutional animal care regulations of Beijing University of Chinese Medicine and were conducted in accordance with AAALAC and IACUC guidelines. Keeping mice in specific pathogen-free conditions in Beijing University of Chinese Medicine. All mice were kept at a controlled room (25 ± 1°C, 45–60% humidity). The allergic asthma model of mice was established by OVA sensitization and atomization inhalation stimulation. Briefly, mice were intraperitoneally injected with 2 mg of OVA (Sigma-Aldrich, Cat#A5503) mixed with 2 mg Imject™Alum Adjuvant (Invitrogen, Cat#77161) and PBS on day 0 and day 14. The challenge phase was from the 21st day to the 25th day after injection, and the mice were atomized with 1% OVA for 30 minutes. Then, the animals were killed, and part of the lung tissue was taken and placed in 4% paraformaldehyde fixed solution for the preparation of histopathological sections. The remaining lung tissues were washed with PBS and stored in - 80°C refrigerator for reverse transcription quantitative real-time polymerase chain reaction (RT - qPCR).

Rabbit anti-PTK6 polyclonal antibody (abs117941), rabbit anti-MAP2K7 polyclonal antibody (abs136429) and rabbit anti-CD46 polyclonal antibody (abs136051) were purchased from Aibixin Biotechnology Co., Ltd. (Shanghai, China). Horseradish peroxidase-labeled goat anti-rabbit IgG antibodies (GB23303) were purchased from Wuhan Servicebio Technology Co., Ltd. After being dewaxed, the lung slices were subjected to antigen recovery with citrate buffer under microwave heating. The slices were cooled down to room temperature and then sealed with 3% bovine serum albumin (BSA) for 30min and were incubated overnight with primary antibody at 4°C. The primary antibodies used were rabbit anti-PTK6 antibody (diluted 1:200), rabbit anti-MAP2K7 antibody (diluted 1:200), and rabbit anti-CD46 antibody (diluted 1:200). The slices were then washed with PBS and incubated with goat anti-rabbit secondary antibody (diluted 1:200) at 37°C for 50 min. After being rinsed with PBS, the slices were visualized with diaminobenzidine and counterstained with hematoxylin.

mRNA was extracted from the lung tissue using a universal RT-PCR Kit (Solarbio Science & Technology Co., Ltd., Shanghai, China) following the manufacturer’s instructions. Samples were treated with DNase and then purified using an RNeasy kit (Qiagen, Hilden, Germany). Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was used as internal reference. PCR primer sequences included the following: ACBD5: forward primer: 5’-TCGCAGGCGAAATTATCTTTG-3’; reverse primer: 5’-GTGCCAACCACTGAGCAATAA-3’, S100A9: forward primer: 5’-CATAAATGACATCATGGAGGACC-3’; reverse primer: 5’-TTGCCATCAGCATCATACACTC-3’, CD46: forward primer: 5’-TGGAAGGCAGTAGCATGGTGAT-3’; reverse primer: 5’-GAGGCTTGGTAGGATGAGTAGGC-3’, MAP2K7: forward primer: 5’-CAATGACTTGGAGAACTTGGGTG-3’; reverse primer: 5’-CGCCGCATTTGCTTAACAG-3’, PTK6: forward primer: 5’-CTGCGTGACTCTGATGAGAAAGC-3’; reverse primer: 5’-AGGTCACGGTGGATGTAATTCTG-3’, GAPDH: forward primer: 5’-CCTCGTCCCGTAGACAAAATG-3’; reverse primer: 5’-TGAGGTCAATGAAGGGGTCGT-3’.

Data were shown as mean ± standard deviation (SD). The differences between the two groups were evaluated by independent sample t-test and nonparametric test. A P < 0.05 was considered statistically significant. Statistical analyses and figures were obtained using IBM SPSS Statistics 23.0 (IBM SPSS Software, NY, USA) and GraphPad Prism Version 8.0 (GraphPad Software, San Diego, CA, USA).

The study involved 531 ARGs. A total of 66 differentially expressed ARGs were observed in mild to moderate and severe asthma. Among these ARGs, SERPINB10 and ATG9B had the highest degree of difference (P < 0.001, Figure 1A ). These ARGs had a PPI network ( Figure S1A ), and the MCODE analysis showed that this network had a key module ( Figure S1B ). Co-expression relationships were explored for 21 differentially expressed ARGs at P < 0.01, and a close correlation between genes was found ( Figure 1B ). These findings suggested that ARGs were involved in the disease process of asthma and might be related to the evolution of clinical symptoms or course of disease.

To investigate the role of ARGs in asthma, we performed the unsupervised consensus cluster analysis on asthma samples based on the expression of 66 ARGs in combination with clinical data ( Figures 1C and 2A, B ). Differences in the expression of 14 ARGs were observed in the two different asthma subtypes (C1 and C2), and the expression of ARGs was downregulated in C1-subtype samples, which were the low-autophagy subtype. The analysis of relevant clinical features showed a significant difference in the use of ICS between patients with C1 and C2 subtypes (P < 0.05). The C1 subtype was higher than the C2 subtype, suggesting the presence of glucocorticoid resistance, and patients with the C1 subtype had more severe symptoms than patients with the C2 subtype ( Figure 2C ). Other key clinical indicators also supported this result. The asthma control questionnaire (ACQ) is the most widely used standardized questionnaire for asthma control and can evaluate the control degree of patients stably and effectively (36–38). The FEV1 reversibility is a marker of airway hyper-responsiveness (AHR) and has traditionally been considered positively correlated with the severity of asthma (39, 40). A combined analysis of two key indicators, i.e., ACQ control and FEV1 reversibility, found that patients with the C1 subtype had higher FEV1 reversibility and more severe poor asthma control than those with the C2 subtype ( Figure 2D ). Even at the same level of FEV1 reversibility, patients with the C1 subtype had lower asthma control than those with the C2 subtype (P = 0.011). These results suggested that ARGs were involved in the disease exacerbation process of asthma, associated with clinical outcomes, and had a good classification function.

To explore the biological characteristics of the C1 subtype, we performed GO analysis on the DEGs of C1 subtype. Results showed that the generation of precursor metabolites and energy, hormone metabolic process, and cilium assembly were the main biological processes; ribosomes and mitochondria were the key cell components; and oxidoreductase activity was the main molecular function ( Figure S1C ). The KEGG analysis revealed that neurodegeneration, amino acid metabolism, and energy metabolism were core pathways in this group ( Figure 2E ). GSEA showed the difference of main signal pathways between C1 and C2 subtypes ( Figures 2F, G ), in which the C1 subtype was predominantly enriched in complement and coagulation cascades, hematopoietic cell lineage, and other pathways. We further constructed the WGCNA network and explored biomarkers in C1-subtype samples on the basis of clinical features. Gene modules were detected on the basis of the TOM, and 28 modules were detected ( Figures 3A, B ). Further analysis of the relationship between modules and clinical features showed that MEdarkslateblue modules had the highest correlation with asthma control (COR = 0.58, P = 0.02). Therefore, the MEdarkslateblue module was selected for subsequent analysis. GO analysis showed that the characteristic genes in the module were involved in the Wnt signaling pathway, positive regulation of nitric oxide biosynthetic process, and extracellular matrix organization, which were closely related to AHR and airway remodeling in asthma ( Figure 3C ) (41). Genes within the same cluster often shared common transcription factors, and we predicted and analyzed the transcription factors of genes in the MEdarkslateblue module and visualized the mutual regulatory relationships between the top 10 transcription factors in Mean Rank ( Figure 3D ). The association of these transcription factors with asthma was demonstrated. For example, histone deacetylase 4 can mediate KLF5 deacetylation to upregulate CXCL12, leading to airway remodeling and promoting the progression of asthma (42).

Based on the DEGs between C1 and C2 subtypes, we identified 12 genes as diagnostic markers for asthma by using the LASSO logistic regression algorithm ( Figure 4A ). Two genes were identified as diagnostic markers for asthma by using the SVM-RFE algorithm ( Figure 4B ). The genetic markers obtained by the two algorithms overlapped to obtain SERPINB10 as the diagnostic gene. The expression of SERPINB10, an ARG, in the C1 subtype was significantly lower than that in the C2 subtype (P < 0.01, Figure 4C ). To assess the immune landscapes of C1 and C2 subtypes, we quantified the level of immune cell infiltration, and the ratios and differences of 22 types of immune cell infiltration between the two subtypes were presented as heatmaps and histograms ( Figures 4D and S1D ). The violin diagram of differences in immune cell infiltration showed that compared with the C2 subtype, the C1 subtype had less CD8⁺ T cell, resting dendritic cell, and resting mast cell infiltrations and more activated mast cell and neutrophil infiltrations (P ≤ 0.05, Figure 4E ). A strong correlation was also observed between several types of immune cells closely associated with asthma. CD8⁺ T cells was significantly negatively correlated with activated mast cells and M0 macrophages, whereas activated mast cells were significantly positively correlated with eosinophils and M2 macrophages ( Figure 5A ). The correlation analysis between SERPINB10 and immune cells showed that SERPINB10 was positively correlated with CD8⁺ T cells (R = 0.34, P = 0.036), T follicular helper cells (R = 0.35, P = 0.03), resting dendritic cells (R = 0.5, P = 0.0013), resting mast cells (R = 0.53, P = 7e−0) and negatively correlated with activated mast cells (R = −0.36, P = 0.028) and neutrophils (R = −0.46, P = 0.0035; Figure 5B ). Finally, to explore the gene expression regulatory network of SERPINB10, a key ARG, we explored the ceRNA mechanism of SERPINB10 and constructed a ceRNA network including 38 lncRNAs and 2 miRNAs ( Figure 5C ). The ERPINB10 gene is regulated by miR-876-5p and miR-2681-3p, and the downregulation of miR-876-5p and miR-2681-3p may lead to the upregulation of SERPINB10 gene in asthma.

We analyzed ARGs and clinical data and found 61 ARGs associated with asthma control ( Figure 6A ). The LASSO Cox regression analysis was performed on ARGs, and a prediction model consisting of five ARGs (i.e., ACBD5, S100A9, CD46, MAP2K7, and PTK6) was further obtained ( Figure 6B ). The risk score for each patient was calculated as follows: Risk score = (−0.539 × ACBD5 expression level) + (0.016 × S100A9 expression level) + (−0.152 × CD46 expression level) + (0.934 × MAP2K7 expression level) + (0.182 × PTK6 expression level). Patients were divided into high- and low-risk groups on the basis of the median cutoff score. PCA based on genomic expression data showed significant distribution differences between the two groups ( Figure S1E ). We found that with increased FEV1 reversibility, the level of asthma control decreased gradually. In addition, poor asthma control was more common in the high-risk group than in the low-risk group. At the same level of FEV1 reversibility, the high-risk group was more likely to be poorly controlled than the low-risk group ( Figure 6C ). This result suggested that risk models based on ARGs could be used as a promising indicator for assessing the control status of patients with asthma. To understand the differences of immune microenvironment characteristics among different risk groups mediated by ARGs, we used ssGSEA to analyze their immune function and immune cell abundance. The high-risk group had higher T cell co-inhibition and lower degrees of immune responses, such as types I and II IFN responses, than the low-risk group, but no statistical difference was observed between the two groups ( Figure 6D ). The immune cell infiltration analysis showed that the high-risk group had a higher degree of B-cell, pDC, and Tfh infiltrations than the low-risk group ( Figure 6E ). These cells were closely related to the occurrence and development of asthma, and the difference was statistically significant (P < 0.05). These findings confirmed the risk-predictive value of predictive models from an immunological perspective. The results of animal experiments were the same as those of the prediction model. ACBD5 and CD46 mRNA were down-regulated in asthmatic mice, while S100A9, MAP2K7 and PTK6 mRNA were up-regulated in asthmatic mice ( Figure 7A ). There was a statistically significant difference in CD46, MAP2K7 and PTK6 (P < 0.05). Further immunohistochemical staining showed that compared with normal mice, CD46 was low expressed in the airway epithelium of asthmatic mice, while MAP2K7 and PTK6 were highly expressed ( Figure 7B ). The analysis of the ceRNA mechanism of the five risk genes in the model showed that S100A9, CD46, miR-328-3p, miR-20a-5p, NEAT1, MALAT1, and XIST occupied an important position and were closely associated with autophagy ( Figure 7C ). For example, miR-20a-5p, which is associated with CD46, is involved in autophagy-induced apoptosis and inflammation in asthma (43), but long-chain noncoding RNA can be used as “molecular sponge” of miR-20a-5p to participate in the regulation of target gene expression (44). Evidence showed that MALAT1 and XIST are involved in the process of autophagy of immune cells (45). Finally, eight potentially important small molecules targeting asthma were screened by the cMAP database, and seven of them had significant correlation with TBXA2R. These molecules were carbacylin, dinoprostone, fluprostenol, GR-32191, iloprost, picotamide, and U-46619, which were considered as PPAR, prostanoid, and thromboxane receptor agonists ( Figure 7D ). Carbacyclin, dinoprostone, fluprostenol, and iloprost were molecularly docked with TBXA2R to confirm their binding ability. Results showed that the affinity between them was all less than −5.0 kcal/mol. Their molecular docking patterns are shown in Figures 7E–H . Small molecules with high affinity might reverse or induce the biological state encoded by specific gene expression markers, thereby playing a potential therapeutic role in asthma.