Infiltrating immunocytes landscape in periodontitis
To explore the differences of infiltrating immunocytes between healthy and periodontitis tissues, CIBERSORT algorithm was conducted to reveal the distribution of 22 immunocyte fractions and we found no matter in healthy or periodontitis tissue the plasma cells took up the largest fraction (Figure 1A). We compared the fractions of 22 immunocytes between healthy and periodontitis, and 13 immunocytes were altered in periodontitis tissue including 5 increasing immunocytes fractions (plasma cell, neutrophils, T cells CD4 memory activated, dendritic cells activated and NK cells resting) and 8 decreasing immunocytes fractions (Tregs, dendritic cells resting, eosinophils, mast cells resting, macrophages M1, B cells memory, T cells follicular helper and T cells CD4 memory resting) in periodontitis tissue (Figure 1B, Table S1). The most downregulated and upregulated immunocytes are Tregs and neutrophils, respectively (Figure 1C). Then, to explore if the variation in immunocytes proportions could be an intrinsic feature that may characterize the individual diversities, we performed PCA analysis and found that the proportions of immune cells displayed distinct clustering and individual differences (Figure 1D). Among them, plasma cells, B cell memory and mast cells resting contributed the most in the healthy and periodontitis tissue variation. Correlation analysis was performed in periodontitis samples and all samples. We found in both periodontitis and all samples that B cell memory and plasma cells were the most negatively correlated , to plasma cells and macrophages M1 (Figure 1E). Further, the differences of infiltrating immunocytes in different clinical subgroups were also investigated including periodontitis subtype, gender, and age by using similar analysis strategies (Figure S1, Table S2-).
To explore and identify the robustness and contributing weight of immunocytes that are related to periodontitis independently, we conducted a series of regression algorithms to adjust the immunocyte interaction effect including univariate logistic regression, LASSO regression, and multivariate logistic regression. 12 periodontitis-related immunocytes were first identified by univariate analysis (Figure 2A, Table S2). Then, LASSO regression was performed to make feature selection and dimension reduction so that the unimportant immunocytes could be excluded. Therefore, we screen out 11 important periodontitis-related immunocytes through LASSO regression (Figure S2). Last, multivariate logistic regression was conducted on the 11 immunocytes evaluate their independency. Only the dendritic cells resting, dendritic cells activated and T cells follicular helper were independently related to periodontitis (Figure 2B).
Identification of immune-related periodontitis subtypes and their characteristics
To identify immune-related periodontitis subtypes, 1214 immune-related genes were compiled, and 702 of them were significantly related to periodontitis (FDR < 0.05). Using the 702-gene panel, 2 distinct immune-related periodontitis subtypes were identified by performing consensus clustering analysis (Figure 3A, Figure S3A-C). Then, 116 genes were selected from all 702 immune-related genes and the sparse structure of the metagenes appeared in the localized patterns of the gene contributions, in which the 116 genes were defined as hub genes for each subtype representing the subtype characteristics (Figure 3B, Figure S3D). There were 45 immune-related genes representing subtype 1 (cluster 1) and 71 genes for subtype 2 (cluster 2). To get a further understanding of the regulation effects among these genes, the PPI networks of the 116 hub genes revealed their relationships and importance in the network., in which CXCR4 of cluster-1 has the most links with other genes and PIK3RI in cluster-2 has the most (Figure 3C, Table S3).
To understand the characteristics and differences between the two immune-related subtypes, GO enrichment analysis and canonical pathway enrichment analysis was performed. Biological processes of GO were conducted and to simplify the results, we grouped categories according to the functional theme, which represents the biological characteristics of the cluster (Figure 4A). The characteristics of cluster-1 were mainly concerning with FC-receptor and immune-associated receptor-mediated immune regulations such as FC/FC-gamma receptor signaling pathway and cell surface receptors. Another feature of cluster-1 is involved in B cell immune processes regulations such as B cell differentiation and B cell activation. As for cluster-2, the IL-6 cytokine feature is predominant such as regulation of IL-production. Then, canonical pathway analysis revealed the significantly enriched or activated pathways in cluster-1 and cluster-2 (Figure 4B-C), which were in consistent with the results above, and we found some famous pathways were shown such as JAK-STAT signaling in cluster-1 and PK3K-AKT in cluster-2.
Associations between immune subtypes and immunocytes
To investigate the diversities between the two immune subtypes from cell level, the 22 infiltrating immunocyte fractions were compared and clinical features were also compared. From the results, the two immune subtypes demonstrated distinct profiles in infiltrating immunocytes, in which 13 immunocytes were different (Figure 5A-B). Distinct infiltrating immunocyte abundances were found between the two clusters of immune periodontitis subtypes. Cluster-1 has more infiltration of activated dendritic cells, plasma cells and neutrophils compared to cluster-2. While for cluster-2, more B cell memory, T cells CD4 memory resting, T cells follicular helper, NK cells activated, monocytes, macrophages M1, macrophages M2, dendritic cells resting, mast cells activated and mast cells resting were upregulated compared with cluster-1. Besides, we noticed cluster-2 was significantly associated with chronic periodontitis and cluster-1 associated with aggressive periodontitis (Figure 5C), while there is no different distribution between two clusters concerning gender and age.
Differentially expressed genes and gene-sets between healthy and periodontitis
To understand the molecular mechanisms by which genes or gene-sets drive the biological reactions in the periodontitis, DEGs analysis and gene-set enrichment analysis were conducted. Using the criterion of FDR < 0.01 and |logFC| > 1, 171 DEGs were identified and 134 of them were up-regulated in periodontitis while 37 of them were down-regulated. These up-regulated genes may be involved in the progress of periodontitis, and among them, the MZB1, TNFRSF17, and IGLL5 are the top three up-regulated genes. The interaction relationships among these genes were also investigated and were presented as a PPI network (Figure 6A). According to outward traversal from a locally dense seed protein and vertex weighting by local neighborhood density, these DEGs were divided into several molecular complexes and were presented by different shapes. Besides, the top ten hub genes among these genes were also identified with the highest links to nodes, indicating they are in the core position of the network and played an essential role in the regulation of the network.
Not only were the individual gene was investigated, but also gene-sets. We mainly focused on four aspects of gene-sets including immunologic signatures, hallmark gene sets, GO biological processes and canonical pathways. Results showed many enriched gene-sets in periodontitis were associated with immune processes or regulations (Figure 6B-E) such as CD4 T CELL VS B CELL DOWN in immunologic signature, IL6 JAK STAT6 SIGNALING in hallmark, POSITIVE REGULATION OF B CELL DIFFERENTIATION in GO and T HEPLER PATHWAY in the canonical pathway. All these four aspects implied immune-related response is essential in periodontitis.
Identification of infiltrating immunocytes related hallmark pathways
To further understand the molecular mechanisms by which hallmark pathways are involved in immunocyte content, we examined the correlation between the fractions of 22 immunocytes and the activity of 50 hallmark pathways. In total, 463 statistically significant correlations between 22 immunocytes and 50 hallmark pathways were observed. Both T cells CD4 naive and B cell memory had 24 negatively correlated signaling pathways, which is the most compared with other immunocytes. While Plasma cells had 21 positively correlated signaling pathways which is the most among all immunocytes. (Figure 7A, Table S4). The links with high correlation coefficients were presented (Figure 7B) such as dendritic cells resting is negatively correlated with hallmark coagulation with -0.56 correlation and plasma cells are positively correlated with hallmark unfolded protein response with 0.57 correlation.
Infiltrating immunocytes related gene module identified by WGCNA
To understand the molecular mechanisms of alteration in infiltrating immunocytes in detail, the genes or gene modules regulating immunocytes or affected by immunocytes need to be uncovered. Therefore, WGCNA, an advanced bioinformatic algorithm, was conducted to find clusters (modules) of highly correlated genes related to external sample characteristics (immunocytes fractions). It is usually used to find a cluster of genes related to some features. 13 was used as network topology for thresholding power considering its relatively mean connectivity power and balanced scale independence (Figure S4A-B). Clustering of module eigengenes was performed in the top 25% variance 4329 genes and 27 modules were identified (Figure 8A-B). Each module was constituted by those highly correlated genes, which have similar functions. Then, the relationship between immunocyte fractions and clinical features were calculated and visualized according to their significance in a heatmap (Figure 8C, Figure S5). In another word, the gene modules related to an immunocyte were identified. For example, the turquoise module was positively correlated with plasma cells the most, as we can see from the heat map and scatter plot (Figure 9A), This suggested the genes in the turquoise module is involved in the regulation of plasma cell population or affected by increased plasma cell fractions. And in another way, it means plasma cell-related genes were identified. The genes with high gene significance and module membership were the key genes in the turquoise module and their interaction relationships were presented in the PPI network so that we could locate the central genes according to links degree (Figure S4C). From the network, we noticed BTK has the most links in the network and is considered a hub gene. GO analysis was also conducted to assess the biological functions of the turquoise module, and it was found these genes were involved in B cell regulation (Figure 9B). This also indicated that our results were reliable from another perspective. Furthermore, a lot of new and unreported plasma cell-related genes were also found (Figure 9C-E, Figure S4D).
An easy-to-use and user-friendly web-tool for periodontitis related immunity analysis
To help periodontitis researchers used the results described and to establish broad access to results in this work as well as other aspects of interest in periodontitis, we established an easy-to-use and user-friendly interactive web-tool for researchers working on periodontitis (http://118.24.100.193:3838/tool-PIA/). We named this platform as Periodontitis Immune Atlas (PIA). PIA provides 6 modules for users: “Genes”, “Pathways”, “WGCNA”, “Correlations”, “Immunocytes” and “Functions”, and all the results in these six modules can be downloaded for usage. The “Genes” module (Figure 10A) provides the DEGs between healthy and periodontitis in several forms. Users can browse the DEG list in a table and sort them by the way they like. Besides, a gene list of immune-related genes with their category is also presented. The specific gene of interest can be visualized by volcano-plot or box-plot. Two genes’ correlation analysis are also provided. The “Pathways” module (Figure 10B) allows users to inquire the activity score of specific pathways and presented them by box-plot. Three commonly used pathway sets were included ( hallmarks, KEGG and GO-BP). The correlation between the two pathways can also be investigated. The “WGCNA” module (Figure 10C) provides analysis for phenotype and gene module. Users can find gene modules related to interesting phenotypes and the gene modules function can also be investigated. In the “Correlation” module (Figure 10D), users can fully explore the relationships among genes, pathways, and immunocytes. “Immunocytes” module (Figure 10E) allows users to visualize the status of immunocytes in different clinical statuses and the correlation between two immunocytes is also included. The “Functions” module (Figure 10F) provides GSEA analysis between different clinical status. Besides, single-gene GSEA can be used as a function predictor for a specific gene in periodontitis. All figures generated in this platform can be downloaded for further analysis and the query results are demonstrated in a comfortable form. The functions provided in the web-tool, which will be continuously updated, may serve as a guide for researchers interested in periodontitis related immunity