Deep Learning Reveals Key Immunosuppression Genes and Distinct Immunotypes in Periodontitis

Background Periodontitis is a chronic immuno-inflammatory disease characterized by inflammatory destruction of tooth-supporting tissues. Its pathogenesis involves a dysregulated local host immune response that is ineffective in combating microbial challenges. An integrated investigation of genes involved in mediating immune response suppression in periodontitis, based on multiple studies, can reveal genes pivotal to periodontitis pathogenesis. Here, we aimed to apply a deep learning (DL)-based autoencoder (AE) for predicting immunosuppression genes involved in periodontitis by integrating multiples omics datasets. Methods Two periodontitis-related GEO transcriptomic datasets (GSE16134 and GSE10334) and immunosuppression genes identified from DisGeNET and HisgAtlas were included. Immunosuppression genes related to periodontitis in GSE16134 were used as input to build an AE, to identify the top disease-representative immunosuppression gene features. Using K-means clustering and ANOVA, immune subtype labels were assigned to disease samples and a support vector machine (SVM) classifier was constructed. This classifier was applied to a validation set (Immunosuppression genes related to periodontitis in GSE10334) for predicting sample labels, evaluating the accuracy of the AE. In addition, differentially expressed genes (DEGs), signaling pathways, and transcription factors (TFs) involved in immunosuppression and periodontitis were determined with an array of bioinformatics analysis. Shared DEGs common to DEGs differentiating periodontitis from controls and those differentiating the immune subtypes were considered as the key immunosuppression genes in periodontitis. Results We produced representative molecular features and identified two immune subtypes in periodontitis using an AE. Two subtypes were also predicted in the validation set with the SVM classifier. Three “master” immunosuppression genes, PECAM1, FCGR3A, and FOS were identified as candidates pivotal to immunosuppressive mechanisms in periodontitis. Six transcription factors, NFKB1, FOS, JUN, HIF1A, STAT5B, and STAT4, were identified as central to the TFs-DEGs interaction network. The two immune subtypes were distinct in terms of their regulating pathways. Conclusion This study applied a DL-based AE for the first time to identify immune subtypes of periodontitis and pivotal immunosuppression genes that discriminated periodontitis from the healthy. Key signaling pathways and TF-target DEGs that putatively mediate immune suppression in periodontitis were identified. PECAM1, FCGR3A, and FOS emerged as high-value biomarkers and candidate therapeutic targets for periodontitis.

Background: Periodontitis is a chronic immuno-inflammatory disease characterized by inflammatory destruction of tooth-supporting tissues. Its pathogenesis involves a dysregulated local host immune response that is ineffective in combating microbial challenges. An integrated investigation of genes involved in mediating immune response suppression in periodontitis, based on multiple studies, can reveal genes pivotal to periodontitis pathogenesis. Here, we aimed to apply a deep learning (DL)-based autoencoder (AE) for predicting immunosuppression genes involved in periodontitis by integrating multiples omics datasets.
Methods: Two periodontitis-related GEO transcriptomic datasets (GSE16134 and GSE10334) and immunosuppression genes identified from DisGeNET and HisgAtlas were included. Immunosuppression genes related to periodontitis in GSE16134 were used as input to build an AE, to identify the top disease-representative immunosuppression gene features. Using K-means clustering and ANOVA, immune subtype labels were assigned to disease samples and a support vector machine (SVM) classifier was constructed. This classifier was applied to a validation set (Immunosuppression genes related to periodontitis in GSE10334) for predicting sample labels, evaluating the accuracy of the AE. In addition, differentially expressed genes (DEGs), signaling pathways, and transcription factors (TFs) involved in immunosuppression and periodontitis were determined with an array of bioinformatics analysis. Shared DEGs common to DEGs differentiating periodontitis from controls and those differentiating the immune subtypes were considered as the key immunosuppression genes in periodontitis.

Results:
We produced representative molecular features and identified two immune subtypes in periodontitis using an AE. Two subtypes were also predicted in the validation

INTRODUCTION
Periodontitis involves the inflammatory destruction of the supporting tissues of teeth. It involves a perturbed local host immune response that is ineffective in countering plaque biofilm microbiota (Meyle and Chapple, 2015). Innate and adaptive immunity work in tandem to counter the infectious challenge posed by oral microbiota, limit the spread of infection, and reestablish periodontal tissue homeostasis (Cekici et al., 2014). This delicately orchestrated process involves the actions of several immune regulatory cell types, including oral epithelial cells (Dutzan et al., 2016), neutrophils (Scott and Krauss, 2011), macrophages, dendritic cells , B cells, and T cells (Gemmell et al., 2002). Regulatory T cells (Tregs) have particularly attracted much recent attention as they engender multiple suppressive mechanisms to inhibit various cells involved in innate and adaptive immunity. The role of Tregs in controlling periodontitis due to their immune-suppressive capabilities has been noted (Alvarez et al., 2018). Immune suppression demands the tandem action of multiple immunosuppression genes, several of which have been demonstrated in the context of periodontal pathology. These include programmed cell death 1 (PD1), PD-Ligand 1 (PD-L1) (Bailly, 2020), and Cytotoxic T-Lymphocyte Antigen4 (CTLA4) (Aoyagi et al., 2000), that function as immune checkpoint inhibitors to modulate B-cells, CD8+ T-cells, and CD4+ T-cells, which can amplify infection and promote tissue damage. Therefore, an immune checkpoint blockade has been proposed as a modality to manage periodontitis. However, existing reports have documented very few immunosuppression genes in the context of periodontitis. It is also recognized that immunosuppressive agents impose a risk for periodontal diseases, inducing gingival overgrowth or other alterations in periodontal tissues (Cota et al., 2010). Immunosuppressive medications for immunerelated disorders such as rheumatoid arthritis or solid organ transplantation are associated with periodontal disease. However, the underlying molecular mechanisms remain unclear, and few genes have been implicated. For instance, specific Human Leukocyte Antigen (HLA)-DR1 genotype is documented to protect from gingival overgrowth induced by cyclosporine A (Cebeci et al., 1996). A more expansive understanding of immune suppression genes that are relevant to periodontal disease pathology can lead the identification of candidate genes and molecular pathways of significant potential translational value. Such data may enable the development of gene and targeted drug therapy for multiple periodontal diseases.
Experimental studies are limited by scale, incomplete or inaccurate existing databases, and the cost-intensive nature of molecular experiments, so approaches that can predict previously unidentified gene functions, enable gene function discovery, and automate the identification of inaccuracies can be very valuable (Chicco et al., 2014). Deep learning (DL) computational frameworks are capable of these. In this regard, an autoencoder (AE), is essentially a dimensionality reduction tool, as the "building block" of DL, comprises of a three-layered unsupervised artificial neural network that performs extraction of representative features (Lee et al., 2009;Wang et al., 2016). The AE has been implemented as a DL framework to predict survival in liver cancer (Chaudhary et al., 2018), breast cancer (Tan et al., 2014), head and neck squamous cell carcinoma (HNSCC) (Zhao et al., 2019), and when applied to RNA-seq data (Xiao et al., 2018) has shown value in generating key features from gene expression data that are linked to clinical outcomes.
To our knowledge, the present study is the first to integrate multi-omics data pertaining to immunosuppression genes in periodontitis using a DL-based AE combined with a support vector machine (SVM) classifier (Ju et al., 2015) confirmed in a validation set, along with an array of bioinformatic analysis, with an aim to identify the most significant immunosuppression genes relevant to the pathogenesis of periodontitis.

Study Design
An overview of the workflow of this study is depicted in Figure 1. In brief, two cohorts of periodontitis datasets (GSE16134 and GSE10334) and immunosuppression genes were included. FIGURE 1 | Overall workflow. The flowchart depicts the autoencoder (AE) architecture and workflow combining deep learning (DL) techniques to identify key immunosuppression genes in periodontitis. Immunosuppression genes related to periodontitis from GSE16134 were applied as input features for an AE. The new transformed features in the bottleneck layer of the AE were clustered into different subtypes using K-mean clustering. Then, based on the clustering labels, we selected the top 100 most related genes from GSE16134 based on ANOVA F values. The input dataset was split at a 60%/40% ratio (training set/test set) to assess the robustness of the AE, using a 5-fold CV. Subsequently, based on the above labels of GSE16134, an SVM classifier was built and further applied for prediction in a validation set (GSE10334). To explore the biological roles of the different identified subtypes, differentially expressed genes (DEGs) and transcription factors (TFs), differential expression analysis, functional enrichment analysis, and construction of TF-target DEGs interaction network were, respectively, applied. Eventually, to identify the immunosuppression genes that might be most pertinent to periodontitis, the overlapping DEGs among the DEGs discriminating disease (periodontitis) and controls and DEGs discriminating the subtypes classified with the AE and SVM models were determined.
First, immunosuppression genes related to periodontitis from GSE16134 were identified and applied as input features to build an AE model. Second, each of the new transformed features in the bottleneck layer of the AE was clustered into different subgroups using K-mean clustering. In addition, based on the clustering labels, we selected the top 100 most related genes from GSE16134 based on ANOVA F values. Data partitioning of the inferring samples of GSE16134 was applied to assess the robustness of the AE, using a 5-fold CV. The samples were randomly split into 5 folds, 3 of which were used as the training set (60%) and the remaining 2 (40%) as the test set. Thereafter, based on the clustering results and the top 100 genes of GSE16134, a SVM classifier was built with a 5-fold CV to identify the optimal hyperparameters, and a validation set (immunosuppression genes related to periodontitis in GSE10334) was applied for SVM to predict the subtypes. To explore the biological roles of the different identified subtypes, differentially expressed genes (DEGs) and transcription factors (TFs), differential expression analysis, functional enrichment analysis, and construction of TF-target DEGs interaction network were, respectively, applied. Finally, to identify the immunosuppression genes that might be most pertinent to periodontitis, the overlapping DEGs among the DEGs discriminating periodontitis and controls and DEGs discriminating the subtypes classified with the DL-based model were determined.

Pre-processing of the Dataset
Transcriptomic data from gingival tissue samples affected with periodontitis and the corresponding controls (GSE16134 and GSE10334) were obtained from the Gene Expression Omnibus (GEO) database of NCBI 1 . Detailed information of the two datasets is listed in Table 1. Immunosuppression genes were obtained from databases DisGeNET 2 and HisgAtlas 3 . From these obtained genes, 1,207 immunosuppression genes related to periodontitis were extracted. Next, the two datasets were stacked, and 1,181 immunosuppression genes' expression profiles were found matching in the two datasets. Subsequently, the two datasets were standardized using the "scale" function in R, setting the parameters as (scale = TRUE and center = FALSE).

Features Transformation
Immunosuppression gene expression profiles of 241 disease samples in the GSE16134 dataset were selected as the input for the AE. The re-coding of the DL algorithms was performed using the Python library "Keras 4 ". An AE is a three-layered neural network consisting of input, hidden, and output layers (Wang et al., 2016), and here an AE with three hidden layers was implemented with 200, 100, and 200 nodes per layer each. One hundred nodes produced by the bottleneck layer were regarded as the new compressed representative features of the data. In accordance with previous research, the AE was set up using the following equations (Chaudhary et al., 2018).
To control overfitting, the penalty values αα and αw (the activity regularizer of layer output) were set to 0.00002 and 0.00001. In addition, the AE was trained using the gradient descent algorithm with 20 epochs and 50% dropout. Here, an epoch is an iteration that indicates the number of passes of the entire training dataset, while the size 20 is one of the appropriate training cycles calculated in the evaluation of the model.

K-Means Clustering to Identify Subtypes of Immunosuppression Genes in Periodontitis
The 100 nodes from the bottleneck-hidden layer were considered as new features for the analysis and were clustered with the K-means algorithm. The optimal number of clusters was determined based on two metrics: Silhouette index (Rousseeuw, 1987) and Calinski-Harabasz index (Calinski and Harabasz, 1974), using scikit-learn package (Pedregosa et al., 2011).

Comparison of AE With PCA Based Clustering
Principal component analysis (PCA), a conventional dimension reduction approach was applied to compare with the AE performance (Chaudhary et al., 2018). The same number (100) of the principal components were set as the features in the 4 https://github.com/fchollet/keras bottleneck layer and clustering performances of AE and PCA were evaluated using the Silhouette index (Rousseeuw, 1987).

Data Partitioning and Robustness Assessment
Data partitioning of the inferring samples of GSE16134 was done to assess the robustness of the model, using a cross-validation (CV)-like procedure, as described in earlier reports (Chaudhary et al., 2018;Zhao et al., 2019). First, the samples were randomly split into 5 folds, 3 of which were used as the training set (60%) and the remaining 2 (40%) as the test set. Using this CV approach,10 new combinations (folds) were obtained. In each, a distinct AE and a classifier were constructed in each training fold and were used for predicting the labels in the test set. Eventually, category labels were inferred using an AE based on all the samples, and these labels were used for predicting labels of the validation dataset.

Supervised Classification
First, the obtained features from GSE16134 were standardized with the "scale" function in R, setting the scale as (center = TRUE and scale = TRUE). Then, the top 100 "most relevant" immunosuppression genes in GSE16134 were selected based on the clustering labels and analysis of variance (ANOVA) F values. Since the top 100 genes were also present in GSE10334 dataset, a complementation test for missing genes was not conducted. Subsequently, based on the labels assigned using GSE16134, a SVM classifier was built and further applied for prediction in a validation set (GSE10334). The "scikit-learn" package (Pedregosa et al., 2011) was used to perform a grid search for the identification of the optimal hyperparameters for the SVM model using a 5-fold CV.

Evaluation of the SVM Classifier
Accuracy and area under the curve (AUC) were selected as two metrics to evaluate the performance of the SVM classifier. The percentage of accuracy was calculated as: Accuracy (%) = Predict number / Test number. A receiver operating characteristic (ROC) curve was plotted for the model using the "pROC" (Robin et al., 2011) and the "ggplot2" packages in R 5 . The AUC is the area under the ROC curve, where an AUC value above 70% is considered acceptable (Mandrekar, 2010).

Differential Expression Analysis
Differential expression analysis was performed for each of the datasets (GSE16134 and GSE10334), to identify genes discriminating between the disease and control samples, using the "Linear Models for Microarray data" ("limma") package in R (Ritchie et al., 2015). Genes with P value < 0.05, and |log FC| ≥ 1 was selected as differentially expressed genes (DEGs). The DEGs with Log FC ≥ 1 was defined as up-regulated DEGs, while the DEGs with log FC ≤ −1 were defined as down-regulated DEGs. Differential expression analysis was also similarly conducted for the classified subtypes. Here, genes with P value < 0.05, and To identify the most critical immunosuppression genes in periodontitis, the DEGs discriminating disease and control samples that overlapped with DEGs discriminating the different subtypes were identified and visualized using a Venn diagram. To evaluate the performance of each such identified gene, a ROC curve was plotted as described earlier.

Functional Enrichment Analysis
The DEGs overlapping in the two datasets (GSE16134 and GSE10334) were identified using the "ClusterProfiler" package in R (Yu et al., 2012). The functions of these DEGs were explored by investigating their enriched Gene Ontology (GO) terms, particularly biological processes (BPs) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. The GO/BP terms and KEGG pathways with P value < 0.05 were regarded as significant functions. The top 30 of the enriched GO/BPs and pathways were chosen to be visualized in a bar plot.
In addition, KEGG pathway analysis was applied to determine the characteristics of different subtypes in GSE16134 and GSE10334 each. KEGG pathways with P value < 0.05 were regarded as significant functions. The top 20 of the enriched pathways were listed and visualized using the heatmap function in R (Galili et al., 2017).

Construction of TF-Target DEGs Interaction Network
TF-target gene interaction pairs were downloaded from multiple databases, including TRRUST 6 , cGRNB 7 , HTRIdb 8 , ORTI 9 , and TRANSFAC 10 . The TFs targeting DEGs overlapping in the two datasets (GSE16134 and GSE10334) were extracted and used for constructing the TFs-target DEGs interaction network. The network was visualized using Cytoscape (Version 3.7.2) (Shannon et al., 2003), and the topological characteristics of the nodes in the TF-target gene network were determined.

Identification of Two Subtypes of Immunosuppression Genes in GSE16134 by AE
The optimal number of clusters was determined based on two metrics: Silhouette index (Figure 2A) and Calinski-Harabasz index ( Figure 2B). The value of the silhouette coefficient is between [−1, 1] and the score near 1 indicates a highly dense clustering. When k = 2, the average silhouette width was nearest to 1 (Figure 2A). Using Calinski-Harabasz index, better performance of clustering depends on a higher score and at k = 2, the score (sum of the squared errors) was the highest (Figure 2B). Therefore, the genes were clustered into two subtypes, defined as S1 and S2.

The AE Performed Better Compared to PCA
The performance of the AE was compared to that of PCA based clustering using Silhouette index. While two optimal clusters were extracted by AE (Figure 2C), six optimal clusters were extracted using PCA (Figure 2D), indicating that the difference between PCA transformed features was minimal, and it was difficult to cluster them effectively. Furthermore, the PCA landing points were concentrated in one zone, and the division was not clear. Therefore, the AE emerged as more effective and accurate in clustering features.

SVM Model and Its Validation
Using a 5-fold CV, the input dataset (immunosuppression genes related to periodontitis from GSE16134) were split at a 60%/40% ratio for the training set and testing set. The SVM model presented an accuracy of 92.78% (Table 2), and the AUC score at 97.72%, above 90% (Figure 2E), supporting the model was efficient in distinguishing between classes and thus reliable in predicting significant immunosuppression genes in the GSE10334 dataset (Mandrekar, 2010).

DEGs Involved in Immunosuppression and Periodontitis
Differential expression analysis was applied to the disease and control samples, as well as the two classified subtypes. A total of 236 DEGs consisting of 48 down-regulated DEGs and 188 up-regulated DEGs were identified from the GSE16134 dataset, while a total of 194 DEGs consisting of 42 downregulated DEGs and 152 up-regulated DEGs were identified   (Table 3). For discriminating the designated subtype labels, a total of 219 DEGs consisting of 85 down-regulated DEGs and 134 up-regulated DEGs were identified in the GSE16134, while a total of 240 DEGs consisting of 95 down-regulated DEGs and 145 upregulated DEGs were identified in the GSE10334 dataset ( Table 4). As shown in the Venn diagram ( Figure 3A), three significant DEGs, Platelet Endothelial Cell Adhesion Molecule Frontiers in Genetics | www.frontiersin.org (PECAM) 1, Fc Gamma Receptor (FCGR) 3A, and FOS were found intersecting and considered as potentially most robust immunosuppression genes related to periodontitis. Each of the three DEGs has an acceptable performance, with an AUC value above 70%, listed in Table 5. The ROC curves of the three genes from GSE16134 and GSE10334 are shown in Figures 3B,C, respectively.

Functional Terms Enriched Among the DEGs
Significantly enriched biological processes and signaling pathways related to the immunosuppressive DEGs were identified from those overlapping between GSE16134 and GSE10334. The immunosuppressive DEGs involved in periodontitis were implicated in biological processes, including T cell activation, regulation of lymphocyte activation, regulation of T cell activation, regulation of cell-cell adhesion, and leukocyte cell-cell adhesion ( Figure 4A). The immune activities were mainly regulated by Th17 cell differentiation, cytokine-cytokine receptor interaction, T cell receptor signaling pathway, Th1 and Th2 cell differentiation, Mitogen-activated Protein Kinase (MAPK) signaling pathway, osteoclast differentiation, and Phosphatidylinositol 3-Kinase (PI3K)-Protein Kinase B (Akt) signaling pathway ( Figure 4B). Most pathways of the two subtypes were evident as distinct in GSE16134 (Figure 5), indicating significant differences between the two subtypes in terms of immunosuppressive activities in periodontitis. This difference was also detected between the two predicted subtypes in GSE10334 (Figure 6). Specifically, subtype S1 of immunosuppressive DEGs in periodontitis from both GSE16134 ( Figure 5A) and GSE10334 (Figure 6A) was mainly enriched in cytokine-cytokine receptor interaction, chemokine signaling pathway, Janus kinase (JAK)-Signal Transducer and Activator of Transcription Protein (STAT) signaling pathway, Hypoxia-inducible Factor (HIF)-1 signaling pathway, and T cell receptor signaling pathway. Of note, subtype S1 from GSE16134 was also enriched in PD-L1 expression and PD-1 checkpoint pathway in cancer ( Figure 5A). Whereas subtype S2 was mainly associated with MAPK signaling pathway, osteoclast differentiation, and infection of virus and E. coli bacteria (Figures 5B, 6B).

Identification of Hub Transcription Factors That Targeted DEGs
The TFs-target DEGs interaction network of the immunosuppression genes in periodontitis is shown in Figure 7, consisting of 197 nodes and 447 edges. Top 30 TFs (Table 6) with the highest degree were considered to represent those most critical to this network. Of these, the top 10 TFs in the network were determined as the hubs, including Androgen Receptor (AR), Hypoxia-inducible Factor (HIF)1A, Signal Transducer and Activator of Transcription Protein (STAT) 5B, and STAT4, which were not only TFs but also up-regulated DEGs, and Nuclear Factor Kappa B Subunit 1 (NFKB1), MYC, JUN, Tumor Protein (TP)53, FOS, and Forkhead Box (FOX) O3, which were not only TFs but also down-regulated DEGs.

DISCUSSION
In this study, we used a DL-based algorithm, the AE, for identifying the pivotal immunosuppression genes relevant to periodontitis. With this approach, we re-constructed multiomics data and produced representative molecular features grouped into two immune subtypes and then built an SVM model based on these, which was confirmed using a validation set. Besides, significant pathways and TF-target DEGs involved in immunosuppression during periodontitis were identified.  Notably, we identified the key characteristics of two immune subtypes of periodontitis. We also identified three "master" immunosuppression genes, PECAM1, FCGR3A, and FOS, as candidate genes central to immune suppressive pathogenic mechanisms in periodontitis.
An AE-based DL approach has demonstrated high efficiency and accuracy in predicting biomarker genes for lung cancer, breast cancer, and HNSCC (Xiao et al., 2018). Akin to these studies, CV results indicated this approach was robust in classifying patients into two subgroups. Furthermore, the AE was more efficient and precise in clustering the distinct features, as compared with the commonly utilized unsupervised ordination method, PCA. In addition, the robustness and reliability of the model were confirmed in a validation set.
The central finding of our study is the identification of three distinct immunosuppression genes, PECAM1, FCGR3A, and FOS, which could be potentially high-value biomarkers or candidate therapeutic targets for periodontitis. PECAM1, also known as CD31, is an immunoglobulin (Ig) gene expressed in various cells, such as endothelial cells (ECs), platelets, and immune cells. PECAM1 is found to be a co-modulator of T-cell immunity (Huang et al., 2017) and a promoter of endothelial junctional integrity (Marelli-Berg et al., 2013). Periodontal pathogens, particularly P. gingivalis, can induce vascular damage through the degradation of PECAM1 (Yun et al., 2005;Farrugia et al., 2020). A protective effect of PECAM1 was also detected in transplant arteriosclerosis (Ensminger et al., 2002). FCGR3A is a member of FCGR families, forming a critical link between humoral and cellular immune responses to periodontal microbiota (Chai et al., 2010;Pavkovic et al., 2018). Previous studies have reported single-nucleotide polymorphisms (SNPs) of FCGR3A (rs396991 and rs4455090) were correlated with periodontitis and might impact susceptibility to periodontitis (Kobayashi et al., 2001;Chai et al., 2010). Besides, FCGR3A polymorphism and the allele rs396991 was identified as an independent susceptibility marker of allograft rejection in patients after organ transplants, highly responsive to natural killer (NK) cells (Paul et al., 2019). FOS was also identified as a significant TF in the study.
Of the top 10 hub TFs, six "leader" immunosuppressive TFtarget DEGs with plausible literature evidence were identified as key to periodontitis pathogenesis and included the downregulated TFs (NFKB1, FOS, and JUN), as well as up-regulated TFs (HIF1A, STAT5B, and STAT4). NFKB1, also termed NF-κB, is a core TF implicated in immune and inflammatory diseases (Tak and Firestein, 2001). Periodontal pathogens can activate NF-κB, and thus inhibition of NF-κB might be a therapeutic target for periodontitis (Ambili et al., 2005). Furthermore, NF-κB is activated in transplanted tissue, and its blockade may be potent in preventing allograft rejection after solid organ transplants, considering the role of NF-κB in T cell activation and differentiation (Molinero and Alegre, 2012). FOS is implicated in periodontitis progression acting via the regulation of T-cell receptor (TCR) signaling (Maekawa et al., 2017). C-Jun (encoded by JUN) signaling is activated by Receptor Activator of Nuclear Factor Kb Ligand (RANKL) and essential for osteoclast differentiation (Ikeda et al., 2004). Activator Protein (AP)-1 is a heterodimer composed of the Fos and Jun subunits, which downregulates osteoprotegerin and is highly expressed in periodontal ligament cells, suggesting their role in bone resorption during periodontitis (Suda et al., 2009). Inhibition of c-Fos/AP-1 by T-5224 (a novel chemical) could attenuate inflammation, T cell proliferation, and allograft rejection in pancreatic islet transplantation (PIT) (Yoshida et al., 2015) and be suggested as a target for immunosuppressive therapy. HIF1A/HIF1, an oxygen-regulated subunit (Corrado and Fontana, 2020), is involved in the immune response of periodontitis, playing a pleitropic role in defending against macrobiotics and facilitating the progression of periodontitis (Wang et al., 2017). HIF1 was also suggested to mediate inflammation and immune responses after organ transplantation, mediating angiogenesis and allograft in the donor organs (Xu et al., 2019). STAT5B and STAT4 are members of the STAT family that play important roles in activating gene transcription through various cytokines. STAT5B and STAT4 can be activated by a variety of cytokines, including Interleukin (IL)12, Type I Interferon (IFNI), IL23, IL2, IL27, and IL35 (Garcia de Aquino et al., 2009;Sanpaolo et al., 2020;Yang et al., 2020), which are prominently involved in mediating immune responses during periodontitis. IFN-γ could stimulate the expression of Indoleamine 2,3-Dioxygenase (IDO)1, a critical immunosuppression protein, in primary human periodontal ligament stem cells (Andrukhov et al., 2017). Thus, evidence suggests STAT5B and STAT4 may mediate immunosuppression during periodontitis. The immunosuppression DEGs in the two subtypes were functionally related to multiple immune-related biological processes and pathways, and the two subtypes were distinct in their regulating pathways. In subtype S1, PD1/PLL1 checkpoint signaling, T cell receptor signaling, and signaling pathways related to immunosuppressive factors, including cytokines, chemokines, Janus Kinase (JAK) -STAT, and HIF1, are found to activate up-regulated TFs, such as HIF1A, STAT4, and STAT5B (de Souza et al., 2012). Whereas the signaling pathways enriched in subtype S2 primarily regulated the MAPK signaling pathway and osteoclast differentiation, as well as the infection of virus and E. coli bacteria, targeting the down-regulated TFs, such as NFKB1, FOS, and JUN (de Souza et al., 2012). Immune response-related pathways were mainly involved in the subtype S1, supporting a hypothesis that periodontitis patients with molecular subtype S1 may be more sensitive to and thus respond comparatively well to the immune-related target therapy.
Considering PDL1/PD1 signaling that characterized the subtype S1, it has been found that peptidoglycans from P. gingivalis can lead to the up-regulation of PDL1 expressed by gingival keratinocytes, as well as the overexpression of PD1 expressed on T lymphocytes (Bailly, 2020). The interaction between PDL1 and PD1 can suppress the initial activation and effector function of T cells and thereby promote the progression of periodontal inflammation . As PDL1-inhibitor has shown significant effects as a cancer therapy in clinical trials (Kim et al., 2020), it may also hold potential as immune therapy for periodontitis patients, especially in the case of immune-compromised patients. The inhibition of the JAK-STAT pathway has been indicated as a potential strategy for immunosuppression therapy, targeting the key cytokines, such as IFNg and IL12 (O'Shea and Plenge, 2012). HIF1A pathway has been found to modulate immunosuppressive molecules, typically VEGF, in periodontitis (Vasconcelos et al., 2016), and tumor microenvironment (El-Sayed Mohammed Youssef et al., 2015). Manipulation of the HIF1A pathway has been proposed as a therapeutic intervention in tumor immunotherapy (Li et al., 2018). The MAPK pathway identified in subtype S2, consists of three family sub-members, extracellular regulated kinases (ERK), c-Jun N-terminal activated kinases (JNK), and p38, and is closely related to osteoblast differentiation (Rodríguez-Carballo et al., 2016). Further, inhibition of p38 may particularly have potential therapeutic value in limiting periodontitis progression at multiple levels of the immune response via its effects on different extracellular stimuli (Kirkwood and Rossa, 2009). Of note, bone resorption, a hallmark of periodontitis, is mainly affected through RANKL, a vital osteoclast differentiation factor (Taubman et al., 2005) and Tumor Necrosis Factor (TNF)-a, majorly activated by MAPK and NF-κB pathways (Ketherin and Sandra, 2018), indicating a key role of these pathways in osteoimmunology.
Altogether, using the DL-based predictive model and bioinformatic analysis, our study provides a predictive and theoretical description of functions and mechanisms relevant to immunosuppression genes active in periodontitis pathogenesis. The validated efficiency and accuracy of the DL-model overcome the bottlenecks of current evidence and suggest new insights valuable for potential translation in therapeutic gene targeting. However, considering our study is the first to apply DL methods in the periodontal disease context, it is expected that further well-designed investigations can validate the model considering other aspects of periodontal disease, where specific and precise associations between clinical parameters and target genes might be identified. One caveat of our study is the lack of phenotype information about the periodontitis cases which were grouped into two distinct immune subtypes. Periodontitis is well recognized as a multifactorial disease, where a disease phenotype may result from multiple factors in a "sufficient cause model" (Heaton and Dietrich, 2012). Distinct "immunotypes" in periodontitis may represent heterogeneity in the core biological mechanisms contributing to disease in different subjects.
A more in-depth understanding of these could support precision medicine approaches in the future. Besides, the possible clinical translation of these results may include multiple directions. For instance, the identification of immunosuppression genes may direct the development of improved topical drugs for delivery at diseased periodontal sites, which could avoid side effects inherent to conventional drugs such as antimicrobials. Also, these findings support a hypothesis that manipulation of the identified immunosuppression genes or selection of the drugs targeting immune checkpoints could be protective against periodontal diseases in patients who have had longtime immunosuppressive therapy, such as those with organ transplantation.

CONCLUSION
The DL-based model applied in this study was reliable and robust in predicting immunosuppression genes in periodontitis. An array of pathways and TF-target DEGs were found to be implicated in the immunosuppressive activity during periodontitis. Three "master" immunosuppression genes, PECAM1, FCGR3A, and FOS, were identified as critical to immune suppression occurring during periodontal pathology. Taken together, the DL model revealed novel insights into the molecular mechanisms underpinning periodontitis and identified key candidate genes for further translation in the context of risk profiling and therapeutic development.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding authors.

AUTHOR CONTRIBUTIONS
WN conceived of the presented idea, designed the overall study workflow, analyzed the data, prepared the figures and tables, authored and reviewed drafts of the manuscript, and approved the final draft. AA prepared the figures and tables, was involved in the proofreading and deep editing, and approved the final draft. ZS analyzed the data, prepared the figures and tables, and was involved in the discussion of the results. AO analyzed the data, prepared the figures and tables, and was involved in proofreading and deep editing. CL, SH, QO, and MZ were involved in the discussion of the results, and also prepared the figures and tables. XL and YD designed the overall study workflow, analyzed the data, and prepared the figures and tables. RH, DZ, and GS devised the main conceptual idea, supervised the whole work, and approved the final draft. GP, YW, and XH supervised the whole project, and approved the final draft. All authors contributed to the article and approved the submitted version.