Autophagy-Related Genes Are Involved in the Progression and Prognosis of Asthma and Regulate the Immune Microenvironment

Abstract

Background: Autophagy has been proven to play an important role in the pathogenesis of asthma and the regulation of the airway epithelial immune microenvironment. However, a systematic analysis of the clinical importance of autophagy-related genes (ARGs) regulating the immune microenvironment in patients with asthma remains lacking.

Methods: Clustering based on the k-means unsupervised clustering method was performed to identify autophagy-related subtypes in asthma. ARG-related diagnostic markers in low-autophagy subtypes were screened, the infiltration of immune cells in the airway epithelium was evaluated by the CIBERSORT, and the correlation between diagnostic markers and infiltrating immune cells was analyzed. On the basis of the expression of ARGs and combined with asthma control, a risk prediction model was established and verified by experiments.

Results: A total of 66 differentially expressed ARGs and 2 subtypes were identified between mild to moderate and severe asthma. Significant differences were observed in asthma control and FEV1 reversibility between the two subtypes, and the low-autophagy subtype was closely associated with severe asthma, energy metabolism, and hormone metabolism. The autophagy gene SERPINB10 was identified as a diagnostic marker and was related to the infiltration of immune cells, such as activated mast cells and neutrophils. Combined with asthma control, a risk prediction model was constructed, the expression of five risk genes was supported by animal experiments, was established for ARGs related to the prediction model.

Conclusion: Autophagy plays a crucial role in the diversity and complexity of the asthma immune microenvironment and has clinical value in treatment response and prognosis.

Keywords: asthma; autophagy-related genes; diagnostic model; immune cell; prognosis.

Figure 1

Expression of autophagy-related genes (ARGs) in asthma. (A) 66 differentially expressed ARGs between mild to moderate and severe asthma. *P < 0.05, **P < 0.01, ***P < 0.001. (B) Co-expression network of 22 ARGs. (C) Consensus clustering cumulative distribution function (CDF) for k = 2–9.

Figure 2

Two different autophagy-related subtypes identified in asthma by unsupervised clustering of 66 ARGs. (A) Relative change in area under CDF curve for k = 2–9. (B) Heatmap of the matrix of co-occurrence proportions for asthma samples. (C) Composite heatmap showing the relationship between the expression characteristics of 14 ARGs and the incidence of nasal polyps, GINA control, ICS dosage, prevalence of allergic rhinitis, gender, and smoking. *P < 0.05. (D) Kaplan–Meier curves of different gene subtypes ACQ Control and FEV1 reversibility. (E) KEGG analysis revealing the key signal pathway of C1 subtype. (F, G) GSEA of C1 and C2 subtypes.

Figure 3

Identification and functional analysis of C1 subtype phenotype-related genes. (A) WGCNA of the C1 subtype to obtain a cluster dendrogram of coexpressed genes. (B) Module–trait relationships for C1 subtypes. Each module contains the corresponding correlation and P-value. (C) GO analysis of genes represented by the MEdarkslateblue module revealing biological processes associated with prognosis. (D) Transcription factors that regulate the expression of genes represented by the MEdarkslateblue module and their interactions.

Figure 4

Screening of autophagy-related diagnostic markers and differences in immune cell infiltration between C1 and C2 subtypes. (A) LASSO logistic regression algorithm to screen diagnostic markers. (B) SVM-RFE algorithm to screen diagnostic markers. (C) Difference of SERPINB10 gene expression between C1 and C2 subtypes. (D) Heatmap of the degree of infiltration of 22 immune cells in C1 subtype samples versus C2 subtype samples. (E) Violin diagram of the proportion of 22 kinds of immune cells. Markers in red indicate significant differences between the two subtypes.

Figure 5

Correlation of key ARGs with immune cell infiltration and its gene expression regulatory network. (A) Heatmap of correlations of 22 immune cells. The size of the colored square represents the strength of correlation, and red and blue colors indicate positive and negative correlations, respectively. A dark color indicates a strong correlation. (B) Mountain diagram showing the correlation between the SERPINB10 gene and resting dendritic cells, resting mast cells, activated mast cells, T follicular helper cells, CD8⁺ T cells, and neutrophils. The lollipop diagram shows the correlation between SERPINB10 gene and 22 kinds of immune cells, and red marks indicate P < 0.05. (C) CeRNA networks involved in the regulation of SERPINB10 gene expression.

Figure 6

The prognostic model constructed by differentially expressed genes of different autophagy patterns can distinguish between high- and low-risk patients with asthma. (A) DEGs associated with asthma control status. Red and green colors represent high- and low-risk genes, respectively. (B) Distribution of LASSO coefficients for DEGs. Tenfold cross-validation for tuning parameter selection in the LASSO regression. Dotted vertical lines are drawn at the optimal values by minimum criteria and 1 − SE criteria. (C) Kaplan–Meier curves of ACQ control and FEV1 reversibility in patients with different risk groups. (D) Differences in the degree of response to each immune function in high- and low-risk groups. (E) Differences in the abundance of infiltrating immune cells in the immune microenvironment between high- and low-risk groups. *P < 0.05, **P < 0.01.

Figure 7

Gene expression regulatory network based on prognostic model and prediction of potential therapeutic drugs for high-risk asthma (N=3, *P < 0.05, **P < 0.01). (A) mRNA expression of 5 risk genes in lung tissues of asthmatic mice and normal mice. (B) The expression of CD46, MAP2K7 and PTK6 in airway of asthmatic mice and normal mice. (C) CeRNA network involved in regulating the expression of five prognostic genes. (D) Correlation between potential therapeutic drugs and corresponding targets. (E) Binding conformation of TBXA2R and carbacyclin (binding energy = −59.29 kcal/mol). (F) Binding conformation of TBXA2R and dinoprostone (binding energy = −34.08 kcal/mol). (G) Binding conformation of TBXA2R and fluprostenol (binding energy = −16.13 kcal/mol). (H) Binding conformation of TBXA2R and iloprost (binding energy = − 64.50 kcal/mol).