Genotype and allele frequency distributions for HLA-DRB1 and FCGR2B were analyzed to identify potential associations with pemphigus vulgaris (PV) susceptibility. The genotype distribution for FCGR2B c.671T > C (I232T) was as follows: in the control group (n = 200), TT = 167 (83.5%), CT = 29 (14.5%), and CC = 4 (2.0%); in the disease group (n = 86), TT = 65 (75.6%), CT = 18 (20.9%), and CC = 3 (3.5%). Hardy–Weinberg equilibrium (HWE) testing in the control group yielded χ² p = 0.0539 and exact p = 0.0695, indicating that FCGR2B genotype distributions were consistent with Hardy–Weinberg equilibrium.

As part of the candidate gene association study (CGAS) framework, single-locus analyses were conducted under dominant, recessive, and additive genetic models for both loci to evaluate potential inheritance patterns and their contributions to PV susceptibility.

Allele-based association analysis revealed multiple HLA-DRB1 alleles significantly associated with PV after Benjamini–Hochberg FDR correction. The 04:02 allele (OR = 30.10, 95% CI = 15.50–61.40, q < 0.001) and 14:01 allele (OR = 5.66, 95% CI = 2.96–11.10, q < 0.001) were identified as the major risk variants, while 16:01 (OR = 0.05, 95% CI = 0.0004–0.34, q = 0.001) and 11:01 (OR = 0.36, 95% CI = 0.14–0.82, q = 0.037) showed protective effects (Table 2).

These findings were consistent across both Fisher’s exact tests and Firth penalized logistic regression models, supporting the robustness of the observed HLA-DRB1 associations.

Visual representations of the HLA-DRB1 association results are provided in Fig. 1 (left-forest plot) and Fig. 1 (right-volcano plot).

Allele-based associations of HLA-DRB1 with PV. Note: The left panel (forest plot) presents odds ratios (ORs) and 95% confidence intervals for individual HLA-DRB1 alleles, derived from Firth penalized logistic regression. The right panel (volcano plot) displays − log₁₀(p) versus log₂(OR), highlighting the magnitude and statistical significance of each allele

The forest plot clearly demonstrates the strong positive effect of the HLA-DRB1*04:02 and HLA-DRB1*14:01 alleles, while the volcano plot highlights their high significance levels and magnitude of risk contribution relative to other alleles. Together, these results underscore HLA-DRB1 as the dominant immunogenetic determinant of PV susceptibility.

In contrast, the FCGR2B c.671 T > C (I232T) variant did not reach statistical significance in either Fisher’s exact or Firth logistic regression analysis under the dominant genetic model (OR = 1.63, 95% CI = 0.83–3.15, q = 0.15). Additive (OR = 1.51, q = 0.15) and recessive (OR = 1.83, q = 0.41) models yielded similar non-significant trends, suggesting a consistent but modest risk direction for the C allele across inheritance patterns. Although the direction of effect suggested a possible risk increase among C-allele carriers, the result did not survive correction for multiple testing. Given its non-significance, FCGR2B results are presented in Appendix A (Table 8 and Fig. 7).

Comprehensive two-locus genotype combination analyses revealed distinct patterns of association between HLA-DRB1 and FCGR2B variants in relation to pemphigus vulgaris (PV) susceptibility. Combination-level findings identified both higher-risk and potentially protective genotype profiles across the two loci. Formal epistasis modeling further clarified the relative contributions of each gene, indicating that their effects were primarily additive rather than non-additive.

Among HLA-DRB1 genotype combinations (diplotypes), several allele combinations showed strong disease associations. The most frequent risk combinations were DRB1*04:02–DRB1*14:01 (frequency = 0.168, p < 0.001, q < 0.01) and DRB1*14:01–DRB1*04:02 diplotypes (frequency = 0.155, p < 0.001, q < 0.01), both conferring markedly elevated risk consistent with single-allele results. Conversely, genotype combinations containing DRB1*16:01 or DRB1*11:01 exhibited protective effects (p = 0.012 and p = 0.035, respectively). The distribution of common HLA-DRB1 genotype combinations (diplotypes) in cases and controls is illustrated in Fig. 2, emphasizing the overrepresentation of 04:02–14:01 among PV cases.

Distribution of common HLA-DRB1 genotype combinations (diplotypes) in pemphigus vulgaris (PV) patients and controls. Note: Bar plot showing the relative frequency of major HLA-DRB1 genotype combinations in PV patients and controls. Genotype combinations containing HLA-DRB1*04:02 and HLA-DRB1*14:01 are more frequent among PV patients (teal), whereas combinations including HLA-DRB1*11:01 or HLA-DRB1*16:01 occur more frequently in controls (red), consistent with their potential protective effects

Two-locus genotype combination analysis integrating HLA-DRB1 and FCGR2B identified several combinations significantly associated with PV. The strongest risk multi-locus combination was FCGR2B–HLA-DRB1 (C–HLA-DRB1*04:02) (frequency = 11.6%, score = 10.04, p = 9.96 × 10⁻²⁴), followed by T–HLA-DRB1*12:01 (p = 4.95 × 10⁻⁸) and additional combinations involving HLA-DRB1*04:02 (p = 1.22 × 10⁻⁷). For clarity, numeric codes used in the analysis were mapped to standard allele nomenclature (see Appendix Table 6). Conversely, specific two-locus genotype combinations involving the FCGR2B T allele and HLA-DRB1*13 variants (HLA-DRB1*13:02, HLA-DRB1*13:03, and HLA-DRB1*13:05) were associated with negative score statistics (p < 0.05), suggesting a protective direction of effect (Fig. 3).

Association of combined FCGR2B–HLA-DRB1 genotype combinations with pemphigus vulgaris. Note: Each bar represents a distinct two-locus genotype combination ordered by −log10(p). The dashed red line indicates the nominal significance threshold (p = 0.05). Taller bars correspond to genotype combinations showing stronger statistical association with disease risk, whereas shorter bars indicate weaker or non-significant associations

These findings indicate a coordinated contribution of antigen presentation (HLA-DRB1) and Fc receptor regulation (FCGR2B) pathways, consistent with a polygenic model of autoimmune susceptibility in PV.

Epistasis modeling (Fig. 4) demonstrated a dominant main effect of HLA-DRB1 risk alleles, with carriers showing a ninefold increase in predicted disease probability (6% → 56%). Although FCGR2B_C_carrier status was associated with a minor additive increase (11–13%), the interaction_term (FCGR2B × HLA-DRB1) was not statistically significant (p = 0.57), indicating additive rather than multiplicative effects. Predicted probabilities closely matched observed disease rates, reflecting good model calibration.

Predicted and observed disease probabilities for the interaction between FCGR2B and HLA-DRB1 risk alleles. Note: Carriers of HLA_risk_alleles (blue bars) show markedly elevated disease probabilities regardless of FCGR2B status, indicating a dominant HLA effect. The FCGR2B_C_carrier variable exerts a small, non-significant additive effect. Open circles represent observed disease rates, closely aligned with model predictions, indicating good calibration and absence of epistatic interaction

To further validate the observed genetic patterns, multivariable Firth logistic regression models were conducted integrating FCGR2B and HLA-DRB1 predictors (Appendix Table 9). The HLA-only model (Model 2) exhibited strong discriminatory ability (AUC = 0.80), confirming the dominant contribution of HLA-DRB1 to PV susceptibility. The combined model including both genes (Model 3) achieved a comparable performance (AUC = 0.81), while inclusion of the interaction_term (FCGR2B × HLA) did not significantly improve fit (p = 0.570, FDR (q) = 0.410).

These regression-based findings further support an additive rather than epistatic relationship between the loci, consistent with the two-locus genotype combination and interaction analyses described above.

To quantify the cumulative impact of HLA-DRB1 and FCGR2B variants on disease susceptibility, both unweighted and weighted Genetic Risk Scores (GRS) were constructed. The unweighted GRS, calculated as the total number of susceptibility alleles carried by each individual, showed a strong association with disease risk (OR = 9.03, 95% CI = 5.33–16.45, p < 0.001). Each additional risk allele substantially increased disease odds, demonstrating a clear cumulative genetic effect. The model exhibited good discriminative ability with an AUC of 0.82 (95% CI: 0.77–0.87), accompanied by a Somers’ D value of 0.64 and a Brier score of 0.15, indicating satisfactory discrimination and calibration (Table 3).

The weighted GRS, which incorporated regression-based β-coefficients to account for variant-specific effect sizes, demonstrated improved predictive performance compared with the unweighted score (OR = 2.74, 95% CI = 2.24–3.44, p < 0.001). The weighted model showed strong discrimination with an AUC of 0.89 (95% CI: 0.85–0.94), accompanied by a Somers’ D value of 0.78 and a Brier score of 0.10, indicating strong discriminative ability and good calibration (Table 3). These findings suggest that incorporating variant-specific effect sizes enhances the predictive capacity of genetic risk modeling.

To further explore potential non-linear genetic relationships, an Explainable Artificial Intelligence (XAI) framework was applied using an XGBoost classifier. The model incorporated four predictors: HLA_risk_allele, FCGR2B_C_carrier, interaction_term, and GRS_unweighted. The optimized model achieved strong predictive performance, with an AUC of 0.88 (95% CI: 0.82–0.95), accuracy of 0.84, sensitivity of 0.95, specificity of 0.79, and an F1-score of 0.75, indicating robust classification performance.

To further assess model stability and potential overfitting, a learning curve analysis was performed by training the model on progressively larger subsets of the training dataset and evaluating performance on an internal validation subset, as shown in Appendix A (Fig. 8). The learning curve analysis indicated that the XGBoost model achieved comparable training and validation AUC values across increasing training sample sizes, suggesting a balanced bias–variance trade-off without clear evidence of overfitting. Although variability was higher at smaller training sizes, performance stabilized as the dataset increased, indicating improved model generalization and robustness.

Model interpretability was assessed using SHapley Additive exPlanations (SHAP) (Fig. 5). The SHAP summary plot identified HLA_risk_allele as the dominant predictor (mean |SHAP| = 1.69), followed by FCGR2B_C_carrier (0.12). The interaction_term showed minimal contribution (0.01), and the contribution of GRS_unweighted was negligible. These results highlight the predominant role of HLA-DRB1 risk alleles in disease classification, with FCGR2B acting as a secondary modulator.

SHAP summary plot showing the relative contribution of genetic predictors to pemphigus vulgaris classification. Note: HLA_risk_allele exhibited the highest mean |SHAP| value (1.69), indicating the dominant influence on model output. FCGR2B_C_carrier showed a smaller but detectable contribution (0.12), whereas the interaction_term had a minimal effect (0.01). The contribution of GRS_unweighted was negligible. Warmer colors indicate higher feature values (greater disease risk), while cooler tones represent lower values

Overall, the explainability analysis confirmed that the strong predictive signal of the machine learning model was primarily driven by the HLA-DRB1 risk allele, consistent with the findings from the regression-based association analyses and two-locus genotype combination results.

Figure 5 further highlights the predominant role of HLA-DRB1 in disease classification, with FCGR2B acting as a secondary genetic modulator. The minimal contribution of the interaction term suggests limited evidence for strong non-linear gene–gene interactions in the XGBoost model. Overall, these findings are consistent with the regression-based association analyses and two-locus genotype combination results, reinforcing the central contribution of HLA-DRB1 risk alleles to disease susceptibility.

Functional annotation and protein–protein association analyses were performed to elucidate the biological significance of the identified genetic associations. The STRING protein association network (Fig. 6A) indicated a statistically supported functional link between HLA-DRB1 and FCGR2B (expected interactions = 0; p = 0.046). This connection reflects shared immune pathways rather than direct physical binding and suggests coordinated involvement of antigen presentation (HLA-DRB1) and Fc-receptor–mediated regulation (FCGR2B). These findings support a complementary functional relationship between the two immune genes, consistent with the additive effects observed in the genetic and epistasis analyses.

Functional annotation and pathway mapping analyses for HLA-DRB1 and FCGR2B. (A) STRING-based protein–protein interaction (PPI) network illustrating a functional (non-physical) association between HLA-DRB1 and FCGR2B. The link reflects integrated evidence from shared immune pathways rather than direct physical binding (expected interactions = 0, p = 0.046). (B) KEGG pathway mapping highlighting immune-related pathways associated with HLA-DRB1 and FCGR2B (e.g., Staphylococcus aureus infection, phagosome, tuberculosis, asthma, autoimmune thyroid disease). These pathways reflect curated functional annotations rather than statistical enrichment. (C) GO Biological Process (BP) annotation emphasizing antigen processing and presentation, Fc-receptor signaling, and regulation of adaptive immune responses. These terms represent ontology-based functional descriptors of HLA-DRB1 and FCGR2B, not enrichment statistics

KEGG pathway mapping (Fig. 6B) highlighted well-established immune pathways associated with HLA-DRB1 and FCGR2B. The mapped pathways included Staphylococcus aureus infection, phagosome, tuberculosis, asthma, type I diabetes mellitus, graft-versus-host disease, and autoimmune thyroid disease. These pathways converge on mechanisms central to antigen uptake, phagocytosis, and MHC class II–mediated signaling, reflecting the known immunological roles of the two genes. As the analysis was based on two genes, these results represent curated pathway annotation rather than statistical enrichment.

GO Biological Process (BP) annotation (Fig. 6C) further underscored biological processes relevant to antigen presentation, Fc-receptor signaling, and adaptive immune activation. Key annotated terms included antigen processing and presentation, regulation of immune effector processes, and positive regulation of adaptive immune responses. These annotations highlight the complementary immune functions of HLA-DRB1 and FCGR2B, consistent with their coordinated roles suggested by the genetic and interaction analyses.

Collectively, these analyses suggest that HLA-DRB1 and FCGR2B contribute to complementary immunogenetic mechanisms integrating antigen recognition and inhibitory FcγRIIb signaling. This interpretation is consistent with the regression, two-locus genotype combination, and XAI-based findings, which support primarily additive genetic effects in PV susceptibility.