Performance exploration of multi-gene panels of alopecia areata susceptibility and drug-binding targets

Abstract

Objective: This study aims to identify potential target genes and therapeutic drugs for alopecia areata (AA).

Methods: Utilizing training and testing data, we evaluated multi-gene panels derived from commonly upregulated genes in publicly available AA patient datasets. The functions of these genes in biological processes were analyzed to identify special multi-gene panels that may play crucial roles in AA. Differences in immune cell infiltration between AA patients and healthy controls were assessed using gene set variation analysis (GSVA) and the Wald test. Signature genes were further validated in specific subsets using single-cell RNA sequence data. Finally, molecular docking and molecular dynamics simulation were conducted to evaluate interactions between protein structures encoded by signature genes and the potential new drug candidates.

Results: When the cut-off value of log₂FoldChage was greater than 1.0, 51 common upregulated genes were identified in the datasets GSE68801 and GSE45512, and the enrichment analysis of biological process indicated the significant involvement of immune cells in AA. The predictive performance of multi-gene panels demonstrated excellent accuracy in pathways related to "regulation of T cell-mediated cytotoxicity" and "cell killing." GSVA and the Wald test demonstrated that the infiltration of T cells and NK cells in AA patients was significantly higher than in healthy controls. Based on single-cell immune cell subsets, we found that within the macrophage migration inhibitory factor signaling pathway, the interactions between NK T cells, CD8 T cells, and melanocytes were observed exclusively in AA patients but not in healthy controls. This indicates that NK T and CD8 T cells may play an important role in the attack on hair follicles via melanocytes. Additionally, we selected several important biomarkers for molecular docking with interacting chemicals, evaluated the stability of drug-protein binding patterns through molecular dynamics simulation, and identified several potential targeted therapeutic agents.

Conclusion: In this study, we screened several key genes associated with immune cells and potential drug-like chemicals that could serve as targeted therapies for AA.

Keywords: alopecia areata; immune cell infiltration; machine learning; module eigen genes; molecular docking; molecular dynamics simulation; multi-gene panel.

FIGURE 1

Gene expression and biological function of signature genes. (A) Volcano plots of differentially expressed genes of datasets GSE68801 and GSE45512. (B) Venn plot of common upregulated genes of datasets GSE68801 and GSE45512. (C) Relative expression heatmap of 51 common genes (AA, alopecia areata; NC, normal control). (D) Top 20 terms in enrichment results of biological function for the 51 common genes. The dot size and color represent the gene count and P value of each term, respectively.

FIGURE 2

Model training and testing results. (A) Prediction model and relative expression heatmap for genes in the enriched term “regulation of T cell-mediated cytotoxicity.” Performance of top four multi-gene panel receiver operative curves. (B) Prediction model and relative expression heatmap for genes in the enriched term “cell killing.” Performance of top four multi-gene panel receiver operative curves. (C) Performance of training and testing cohorts for 97 machine learning models.

FIGURE 3

Differences in immune cell infiltration and biological function of gene modules of interest. (A, B) Significance tests and fold changes for immune cell infiltration [(A) CIBERSORT and (B) TCIA] scores of gene expression sets from the scalp of AA patients and healthy controls in datasets GSE68801, GSE45512, and GSE80342. (C) Biological functions of two gene modules (negatively correlated brown module and positively correlated saddle brown module), which focused on skin cell development and immune cell activation, respectively.

FIGURE 4

Pearson correlation results of module eigen genes and immune cell infiltration scores. (A, B) Heatmap of Pearson correlation results of module eigen genes and immune cell infiltration scores. [(A) CIBERSORT and (B) TCIA; *P < 0.05, **P < 0.01, ***P < 0.001].

FIGURE 5

Immune cell interactions in single cell RNA-seq cell subtypes. (A) Cell type atlas of six alopecia areata patients and two healthy controls in the single-cell RNA sequence dataset GSE233906. (B) Difference in cell type infiltration on the basis of top 50 biomarkers in each cell cluster. (C) Predictive performance of gene-panel with 16 genes in pooled datasets GSE68801, GSE45512, and GSE80342. (D) Comparative network of cell–cell chat (including CD4 T, CD8 T, NK T, mac/mono/DCs, and melanocyte) signaling pathways for macrophage migration inhibitory factors.

FIGURE 6

Details of ligand–residue interactions in drug-binding pockets.

FIGURE 7

Evaluation of molecular docking and drug-binding pocket. (A) Details of ligand–residue interactions in drug-binding pockets. (B) Crystal structure of the optimal molecular docking position in each available pocket. (C) Top three chemicals with the highest affinity score for protein molecular encode by each signature gene.

FIGURE 8

Stability evaluation of molecular dynamics simulation for the binding mode of methotrexate and HLA-DRA protein. (A) RMSD changes in the two binding pocket modes. (B) Three-dimensional heatmap of the cumulative distribution of free energy landscapes (FELs) over RMSD variations in a 300-nanosecond molecular dynamics simulation. (C) RMSF of all amino acid residues in each binding pocket mode. (D) SASA changes in the two binding pocket modes. (E) Binding conformation of methotrexate and HLA-DRA protein complex. (F) Free energy diagram of solvent stability and contribution of solvation to binding free energy. The dark-green area and violet grid represent free energies less than −0.2 kcal/mol/ų and greater than 0.2 kcal/mol/ų, respectively, at a distance of 5Å from the ligand. (G, H) Interactions of methotrexate with binding residues in the crystal structure (G) and schematic diagram (H) of the HLA-DRA protein.