Landscape of the infiltrating immune cellsAfter removing samples with CIBERSORT P > 0.05, 218 samples, including 73 rejected and 145 non-rejected samples in the GSE21374 dataset, were reserved for CIBERSORT analysis. As shown in Fig. 2, the histogram and heatmap clearly show the composition and proportion of infiltrating immune cells in each sample. Moreover, correlation analysis revealed that the proportion of most immune cells was related to other immune cells; for example, the proportion of plasma cells was negatively correlated with the proportions of M0 macrophages, M1 macrophages, M2 macrophages, CD4 memory activated T cells, gamma delta T cells, CD8 T cells, and follicular helper T cells, whereas the proportion of CD8 T cells was positively correlated with the proportion of naive B cells, regulatory T cells (Tregs), CD4 memory activated T cells, gamma delta T cells, and M1 macrophages (Fig. 2B C). Furthermore, differential analysis showed that the proportions of memory B cells, plasma cells, CD8 T cells, activated memory CD4 T cells, follicular helper T cells, gamma delta T cells, M1 macrophages, and resting mast cells were markedly different between rejection and non-rejection samples (Fig. 2D). The proportions of memory B cells, plasma cells, and resting mast cells were downregulated in the rejection samples compared to those in the non-rejection samples, but the proportions of follicular helper T cells, CD8 T cells, M1 macrophages, activated memory CD4 T cells, and gamma delta T cells were upregulated (Fig. 2D). Thus, we speculated that the rejection of kidney transplantation might be mainly related to memory B cells, plasma cells, CD8 T cells, activated memory CD4 T cells, follicular helper T cells, gamma-delta T cells, M1 macrophages, and resting mast cells.Fig. 2Composition and proportion of infiltrating immune cells in rejection and nonrejection samples revealed by CIBERSORT: (A) Distribution of immune cell ratios in the two sample types; (B-C) Heatmap of the correlations between the 22 immune cell types; (D) Analysis of the differences in immune cell infiltration in the two subtypes_box plot. Each p-value is written above the box plots (NS: P > 0.05, *: P ≤ 0.05, **: P ≤ 0.01, ***: P ≤ 0.001, and ****: P ≤ 0.0001).Identification of infiltrating immune cell-related genes based on WGCNAWGCNA was performed using the GSE21374 dataset to screen infiltrating immune cell-related genes. First, sample clustering analysis suggested that all 284 samples in the GSE21374 dataset could be used to construct a weighted gene co-expression network (Fig. 3A). Subsequently, when β was set as 5 (R2 = 0.85) and modules with a dissimilarity of module eigengenes less than 0.3 were merged, a weighted gene co-expression network including 13 modules was constructed (Fig. 3B-C). The association between these 13 modules and eight differentially infiltrating immune cells (memory B cells, plasma cells, CD8 T cells, activated memory CD4 T cells, follicular helper T cells, gamma delta T cells, M1 macrophages, and resting mast cells) is shown in Fig. 3D. The purple ME module was associated with multiple immune cells (Fig. 3D and Supplemental Fig. S1). Furthermore, we calculated the correlation between genes in the ME purple module and eight differentially infiltrating immune cells and identified 270 infiltrating immune cell-related genes by setting the correlation coefficient > 0.8 and P < 0.05 (Supplemental Fig. S2).Fig. 3Identification of modules associated with infiltrating immune cell-related using the weighted gene co-expression network analysis (WGCNA). (A) Sample clustering map; (B) Analysis of the scale-free index and mean connectivity for various soft-threshold powers; (C) Dendrogram of all differentially expressed genes clustered based on the measurement of dissimilarity (1-TOM). The color band shows the results obtained from the automatic single-block analysis; (D) Heatmap of the correlation between the module eigengenes and traits of infiltrating immune cell. We selected the ME purple module for subsequent analysis. TOM topological overlap matrix, ME module eigengene.Protein-protein interaction (PPI) network and biological functions of infiltrating immune cell-related genesTo investigate the interactions of 270 infiltrating immune cell-related genes at the protein level, a protein-protein interaction (PPI) network was constructed. As illustrated in Fig. 4A, CD8A, TYROBP, and ITGB2 interacted strongly with other proteins. The results of Gene Ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis are shown in Fig. 4B and C. Notably, these 270 infiltrating immune cell-related genes were mainly involved in immune-related biological processes (BP), cellular components (CC), molecular functions (MF), and Kyoto Encyclopedia of Genes and Genomes pathways. For example, immune-related BPs related to 270 infiltrating immune cell-related genes included T-cell activation, adaptive immune responses, and leukocyte proliferation (Fig. 4B). Moreover, immune-related CCs that were involved in 270 infiltrating immune cell-related genes included the external side of the plasma membrane, mast cell granules, and T cell receptor complexes (Fig. 4B). Immune-related MFs involved in 270 infiltrating immune cell-related genes included interleukin-15 receptor activity, MHC protein complex binding, and T cell receptor binding (Fig. 4B). Furthermore, immune-related KEGG pathways that were involved in 270 infiltrating immune cell-related genes included natural killer cell-mediated cytotoxicity, chemokine signaling pathways, Th1 and Th2 cell differentiation, and B-cell receptor signaling pathways (Fig. 4C). These results suggest that these 270 infiltrating immune cell-related genes may play key roles in the rejection of kidney transplantation by regulating the immune response.Fig. 4PPI network and biological functions of infiltrating immune cell-related genes. (A) Protein-protein interaction (PPI) network (B) Biological process (BP), cellular component (CC), and molecular function (MF) of Gene Ontology (GO) enrichment; (C) Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. The size of the bubble indicates the strength of the P value.Identification of prognostic genesTo further identify genes related to the prognosis of patients with kidney transplantation from 270 infiltrating immune cell-related genes, Kaplan (K‒M) survival analysis was performed on the GSE21374 dataset based on the time information of graft loss for kidney transplantation. Interestingly, we found 212 genes among 270 infiltrating immune cell-related genes that were associated with the prognosis of kidney transplantation patients (Supplemental Table S1).Identification of hub genes by machine learning algorithmsTo further identify hub genes related to rejection of kidney transplantation from the 212 prognostic genes, Extreme Gradient Boosting (XGBoost) and Least Absolute Shrinkage and Selection Operator (LASSO) algorithms were performed on the GSE21374 dataset. First, max_depth = 3, learning_rate = 0.1, gamma = 0, and n_estimators = 100 were identified using 10-fold cross-validation and a grid search in the training set. Recursive feature elimination cross-validation (RFECV) feature selection revealed that the parameter step and Cross Validation (CV) should be set as 1 and 2, respectively, under which a module including 43 genes was identified in the training set (Fig. 5A-B). Moreover, the LASSO algorithm acquired a model including seven genes and showed good efficiency in distinguishing rejection and non-rejection samples in the training and testing sets (Fig. 5C-F). Finally, by overlapping genes obtained by the XGBoost and LASSO algorithms, a total of five hub genes, including WARS, CD8A, CRTAM, GBP2, and VAMP5, were identified (Fig. 5G).Fig. 5Identification of hub genes by machine learning algorithms. (A-B) Extreme Gradient Boosting (XGBoost) regression analysis for screening the characterized gene; (C-D) Least Absolute Shrinkage and Selection Operator (LASSO) regression for variable selection; (E-F). Receiver Operating Characteristic (ROC) curves analysis of the diagnostic model in training and validation dataset; (G) Venn diagrams to obtain overlapping genes for XGBoost and LASSO algorithms.Assessment of the value of hub gene prognosis and diagnosisWe analyzed the relationship between the five hub genes and OS and found that the group with lower expression of these genes had a better prognosis than the group with higher expression of these genes (Fig. 6A). To further analyze the ability of these five hub genes to distinguish rejection samples from non-rejection samples, Receiver Operating Characteristic (ROC) analysis of the GSE21374 and GSE36059 datasets was performed separately. Surprisingly, we found that both Areas Under the Curve (AUC) values for each gene in these two datasets were above 0.7 (Fig. 6B-C), which indicated that CD8A, CRTAM, GBP2, WARS, and VAMP5 might be used as diagnostic biomarkers of kidney transplant rejection.Fig. 6Assessment of the value of hub gene prognosis and diagnosis. (A) Kaplan-Meier survival analysis of the five hub genes; (B) Receiver Operating Characteristic (ROC) curves of hub genes in the training dataset; (C) ROC curves of hub genes in the validation dataset.Correlation among hub genes and differentially infiltrating immune cellsFirst, we observed a correlation among CD8A, CRTAM, GBP2, WARS, and VAMP5 and found that there was a positive correlation between the expression of CD8A, CRTAM, GBP2, WARS, and VAMP5 (Fig. 7A-C). Notably, CD8A, CRTAM, VAMP5, and WARS were positively correlated with CD8 T cells, CD4 memory activated T cells, follicular helper T cells, gamma delta T cells, and M1 macrophages but were negatively correlated with plasma cells, memory B cells, and resting mast cells (Fig. 7D). Moreover, GBP2 was positively correlated with CD8 T cells, activated memory CD4 T cells, follicular helper T cells, gamma delta T cells, and M1 macrophages but negatively correlated with plasma cells and memory B cells (Fig. 7D). Therefore, WARS, CD8A, CRTAM, GBP2, and VAMP5 may affect the rejection of kidney transplants by regulating the infiltrating levels of memory B cells, plasma cells, CD8 T cells, activated memory CD4 T cells, follicular helper T cells, gamma delta T cells, M1 macrophages, and resting mast cells.Fig. 7(A) Correlation analysis of 5 overlapping characterized genes; (B) Correlation analysis of 5 overlapping characterized genes p-value; (C) Correlation analysis of 5 overlapping characterized genes; (D) Correlation of intersecting single genes of characterized genes with differential immune cells. Pearson correlation analysis was used to verify the correlation between hub genes and immune cell infiltration.Verification of the expression levels of hub genesTo verify the expression levels of CD8A, CRTAM, GBP2, WARS, and VAMP5 in the rejection samples, we first observed the expression levels of CD8A, CRTAM, GBP2, WARS, and VAMP5 between the rejection and non-rejection samples in the GSE21374 and GSE36059 datasets. Interestingly, the expression levels of CD8A, CRTAM, GBP2, WARS, and VAMP5 were upregulated in the rejection samples compared to the non-rejection samples in the GSE2137 and GSE36059 datasets (Fig. 8A,B). Subsequently, we observed the expression levels of CD8A, CRTAM, GBP2, WARS, and VAMP5 in the non-rejection, ABMR, and TCMR samples in the GSE36096 dataset. In ABMR and TCMR samples, the expression levels of CD8A, CRTAM, GBP2, WARS, and VAMP5 were significantly upregulated compared with those in non-exclusive samples (Fig. 8C). In addition, we verified the expression levels of CD8A, CRTAM, GBP2, WARS, and VAMP5 in the clinical samples. Consistently, quantitative real-time PCR also revealed that CD8A, CRTAM, GBP2, WARS, and VAMP5 were upregulated in the rejection samples compared to the non-rejection samples (Fig. 8D). These results further illustrate that CD8A, CRTAM, GBP2, WARS, and VAMP5 play key roles in kidney transplant rejection.Fig. 8(A) GSE21374 Characterized Gene Intersection Single Gene Expression Boxplot; (B) GSE36059 Characterized Gene Intersection Single Gene Expression Boxplot; (C) Quantitative real-time fluorescence PCR and differential analysis showed that CD8A, CRTAM, GBP2, WARS and VAMP5 were up-regulated in rejection samples compared to nonrejection samples.The Wilcoxon rank-sum test was employed to assess the discrepancy between different samples. Each p-value is written above the box plots (NS: p > 0.05, *: p ≤ 0.05, **: p ≤ 0.01, ***: p ≤ 0.001, and ****: p ≤ 0.0001).
