A total of 2,611 human microarray datasets from tissue biopsies of kidney allografts were obtained from the Gene Expression Omnibus (GEO) database ( Supplementary Table 1 ). According to the Banff standard, samples diagnosed as “non-rejection” are classified as “stable state (STA)”, and samples diagnosed as “antibody-mediated rejection (ABMR)” and “T-cell mediated rejection (TCMR)” are classified as “Mixed state (Mix)”. All other data diagnosed as “borderline” are excluded. The data used are all samples from patients diagnosed as ABMR, TCMR, acute rejection (AR), stable state (STA), chronic rejection (CR), and mixed state (Mix) (2, 3, 8, 13–21). The single-cell datasets of samples were downloaded from GEO (12, 22) (GSE145927 and GSE109564), which were collected from kidney biopsies with a diagnosis of acute ABMR and acute Mix.

Microarray datasets were re-annotated to unify the gene names corresponding to each probe. We aligned all probe sequences from 2,611 microarray datasets to the FASTA file of hg38 using bowtie2 (23) and annotated them using bedtools (24). The samples from each series record of GEO were standardized using limma (25). Two expression matrices were prepared. The first one was acquired by directly merging the matrices from each GSE based on the intersection of genes without removing the batch effect and normalization. The second matrix was progressed with Combat to remove batch effects. Batch effect was removed according to the GPL number and the company category. All the combined data of the second matrix were normalized by the logarithm of 2.

The raw gene expression matrices of scRNA-seq datasets from all renal samples were merged and converted to an Anndata object using the Python package Scanpy (version 1.4.4) (26). Cells that expressed less than 500 genes and genes detected in less than 3 cells were filtered out. Potential doublet cells were then detected and filtered by applying the Python package scrublet (version 0.2) (27) for each sample. Next, doublet-dominated sub-clusters were checked to ensure a low doublet rate in all populations using the previously described method (28). The gene expression levels were normalized by the total UMI count per cell (1e4) with data being log-transformed. The interferences arising from cell cycling genes were eliminated by using the regress_out function of the Scanpy package. Then, highly variable genes (HVGs) in gene expression matrices were identified for further analysis using the highly_variable_genes function of the Scanpy package. Finally, the batch effect was eliminated using the Python package bbknn (version 1.2.0) (29). The dimensionality of HVGs was primarily reduced by principal component analysis (PCA). The first 40 principal components were further summarized by UMAP for dimensionality reduction using the default setting of the UMAP function of the Scanpy package. Cells were clustered with the Leiden algorithm using the leiden function of the Scanpy package. Cell-specific gene markers across all cell types were identified with the get_DEG_single function of Python package PLOGS (https://github.com/ZhangHongbo-Lab/PLOGS) that we developed, with parameter ratio ≥ 0.5 and q-value ≤ 1e-30.

To reclassify the samples of human microarray datasets, we used Scanpy (version 1.6.0) to read the first expression value data and the log1p function of Scanpy to perform logarithmic calculations. Then, HVGs were calculated based on each batch and screened with parameter batches ≥ 1. The second expression value data was screened using the HVG, which was the third expression matrix. The dimensionality of HVGs was primarily reduced by PCA. The first 40 principal components were further summarized by UMAP (30) dimensionality reduction using the default setting of the umap function of the Scanpy package and clustered with the Leiden algorithm using the leiden function of the Scanpy package. In order to make the data classification results clearer and more credible, the batch effect was eliminated again using the Python package bbknn (version 1.2.0). The cluster-specific gene markers were identified with the get_DEG_single function of the Python package PLOGS with parameter ratio ≥ 0.5 and q-value ≤ 1e-30. The reclassified clusters were annotated by differential expression gene groups of six rejection states.

Scissor (31) was performed to identify the cell subpopulations most highly associated with the states of reclassified clusters in bulk RNA-seq data. All of the states of reclassified clusters were merged, and each cell corresponded to the rejection state with the largest positive correlation coefficient.

The expression matrices of the marker genes of each cell type were regarded as the cell characteristics, and CIBERSORT (32) was performed to analyze the cell ratio of the six states using the 2,611 microarray datasets. The R package pheatmap was used to visualize the results, and the cell type was set to scale for horizontal comparison. Then, we output the numerical matrix after the scale, calculated the power function matrix corresponding to the numerical matrix with 2 as the base, and rounded it, selecting the cell type with the value greater than 1 as the potentially important cell type. This value was used as the proportion of the cell type in the corresponding rejection state. The score_genes function of the Scanpy package was used to show the degree of gene expression in different states.

Co-expression network analysis was performed using weighted correlation network analysis (WGCNA, R package) (33). All genes were selected as input matrix. The co-expression network was constructed by the automatic construction function with the parameter power 10. Co-expression network was visualized by Cytoscape (34). IRegulon (35) was used to predict the transcription factors.

NicheNet (36) was used to predict upstream regulatory networks that drive STAT4. All expressed genes in T cells were used as the background of genes. Genes were considered as expressed when they have nonzero values in at least 10% of the cells in a cell type. Here, all ligands were adopted to determine signaling paths between ligands and STAT4.

Putative STAT4-regulated genes excluding those with insignificant correlations were applied to least absolute shrinkage and selection operator (LASSO) regression to identify critical prognostic genes and construct a diagnostic model with low variance and strong universality. The dataset GSE21374 was divided into training and validation sets by 6:4 randomization without replacement. The package glmnet was used in R version 4.1.0.

In the current clinical diagnosis of kidney transplant prognosis, allograft rejection is usually classified into six progression statuses based on divergent histology manifestations: stable (STA), antibody-mediated (ABMR), acute (AR), chronic (CR), T-cell mediated (TCMR), and mixed TCMR with ABMR (Mix) rejection. To evaluate the correlation between state from current clinical diagnosis and transcriptome and understand the cellular and molecular mechanisms that lead to kidney transplant rejection, we first applied the current clinical criteria and analyzed microarray data of 2,611 kidney biopsies (2, 3, 8, 13–21) from patients receiving renal transplantation ( Supplementary Table 1 ). Through PCA, we found that allograft samples with different clinical diagnoses were mixed and randomly distributed in two separate groups ( Supplementary Figure 1A ), suggesting the discrepancy between traditional clinical diagnosis and transcriptomic heterogeneity of disease status. At the gene expression level, different clinical diagnoses were also difficult to be recognized by differentially expressed genes (DEGs) as compared between each rejection state and the stable condition ( Supplementary Figures 1B, C ). Similar results were also observed in other previous transcriptomic analyses (2, 17, 37). These results indicated that the traditional clinical classifications might not accurately distinguish patients with different rejection states at the molecular level.

To precisely refine the rejection status at the molecular level, we constructed a classification pipeline based on UMAP and Leiden ( Figure 1A ). Unsupervised clustering of 2,611 samples from kidney transplantation yielded 6 main rejection states ( Figure 1B ). As expected, each status was distinguished by its signature gene sets and characterized by specific Gene Ontology (GO) functions ( Figure 1C ; Supplementary Table 2 ). Samples in the stable state (STA) showed significant enrichment in renal homeostasis and detoxification function ( Supplementary Figure 2 ) representing a relevant homeostatic condition after receiving transplantation. In addition, samples in fibrosis state (Fib) indicated allografts suffering from fibrotic diseases. The population of inflammatory state 2 (Infla2) seemed to represent a status suffering from dysfunctional leukocytes. With enrichment of tube morphogenesis, renal insufficiency, and humoral immune function, both the progressive state 1 (Prog1) and 2 (Prog2) were characterized to be the activated progression stages of rejection.

Interestingly, the population of inflammatory state 1 (Infla1) was not only correlated with graft-versus-host disease, but also predominantly enriched in immune activation responses with more than 70% of samples in Infla1 showing apparent rejection phenotypes, representing the severest rejection status ( Figure 1C ; Supplementary Figures 2 , 3 ). To further verify its rejection-triggering feature, we mapped these new classifications and phenotypes on human and mouse renal allografts (37) (https://www.ualberta.ca/medicine/institutes-centres-groups/atagc/research/gene-lists, Supplementary Table 3 ). Transcript sets that were reduced after injury or rejection (termed “reduced after injury” and “reduced >90% in day21 mouse allografts”) were enriched in all rejection states that specifically excluded Infla1, suggesting that these states, to some extent, still maintained stability ( Figures 1D, E ). However, the injury- and repair-induced transcripts were highly enriched in Infla1 and, to a lesser extent, Fib, indicating that a part of allograft samples in Fib and Infla1 were damaged ( Figure 1F ).

Previous studies showed that immunoglobins, effective CD8⁺ T cells, and cytotoxic molecules such as GZMB and IFN-γ increase the risk of graft failure (38–42). To identify the group with a high risk of graft failure, we selected transcripts that not only were used in clinical diagnosis but also represented a high risk of allograft loss, and calculated the overall expression scores ( Supplementary Table 3 ). Interestingly, all four HR transcript sets showed the highest expression scores in Infla1, further indicating that Infla1 was the HR status prone to allograft loss ( Figure 1G ). Considering all of the above lines of evidence, we identified Infla1 as the HR state.

To uncover key cell types with a significant impact on HR state, we first collected scRNA-seq datasets (12, 22) to analyze all cell types present in kidney rejection samples. Unsupervised clustering of the scRNA-seq data from 3 patients identified 11 main cell types defined by signature genes ( Figure 2A , left panel; Supplementary Figures 4A, B ; Supplementary Table 4 ). Using the Scissor (31) toolkit, we assigned all cells with each of the rejection states ( Supplementary Figure 4A ) and further applied 23,082 positively relevant cells to better illustrate the relationship between specific cell types and rejection states ( Figure 2A , right panel). Notably, we identified type I (M1) and type II (M2) macrophages and T cells were strikingly aggregated in HR. Previous knowledge recognized that macrophages and T lymphocytes were the dominant cell types infiltrating acutely rejecting grafts (43). T lymphocytes are central in promoting transplantation rejection and organ damage through allorecognition of foreign antigens and effector responses (44) ( Figure 2B ). We further performed CIBERSORT (32) to predict relative ratios of each cell type in rejection states and also revealed the immune-related cell types including macrophages, T cells, and B cells highly aggregated in HR ( Supplementary Figure 4C ). These results revealed the cell-type characteristics of each rejection state, which were strikingly beneficial for clinical diagnoses. Interestingly, the enrichment of T cells in HR had more significant differences, showing that changes in the amount of T cells were much stronger than other immune cell types in transcriptional datasets ( Supplementary Figure 4C ). Indeed, in further distinguishing cell types highly related to HR by applying HR transcript sets to scRNA-seq, we observed that immune cells including macrophages, T cells, and B cells were enriched in the HR state, while T cells appeared to be specifically involved, indicating that T cells were more recruited in driving HR rejection progression ( Figure 2C ).

To explore the potential molecular mechanisms mediated in T cells that lead to HR rejection, WGCNA (33) was performed to seek out the gene co-expression networks ( Supplementary Figures 5A–C ). Eighteen gene modules (labeled with colors, such as MElightcyan) were generated by calculating the correlation between total genes and the allograft samples ( Supplementary Figure 5B ). The number of significant correlation coefficients between gene modules and the newly defined rejection states was far greater than the classification based on traditional clinical diagnoses, suggesting the power of this new classification ( Figure 3A ; Supplementary Figure 5C ). Additionally, the gene module MEblack has the biggest correlation with the HR state ( Figure 3A ). To focus on modules most relevant to HR, we selected MEblack and identified that a large majority of its hub genes were regulated by signal transducer and activator of transcription 4 (STAT4), which was detected by IRegulon (35) ( Figure 3B ; Supplementary Figure 6 ). GO analysis of these hub genes regulated by STAT4 showed highly activated immune functions, including leukocyte activation and regulation of lymphocyte proliferation ( Figure 3C ), indicating that STAT4 mediated high levels of immune responses in HR rejection. Indeed, upon analyzing an external microarray dataset GSE21374 (45) from renal allograft, we observed that patients with a higher expression of STAT4 showed poorer allograft survival in renal transplantation ( Figure 3D ). Intriguingly, STAT4 was not only strikingly expressed in HR ( Figure 3E ), but also significantly enriched in T cells, illustrating that STAT4 mainly conducted HR rejection in T cells ( Figure 3F ). It is well documented that STAT4 is a member of the STAT family, which are identified as the major components of DNA-binding proteins that activate gene transcription in response to a variety of cytokines (46, 47). It contributes to the differentiation and proliferation of both Th1 and Th17 cells, which are also crucial effectors in chronic inflammatory disorders (48). Therefore, highly correlated with the development of autoimmune diseases (47), STAT4 has a large potential to be a key regulator of graft-rejection activation.

We next investigated whether hub genes regulated by STAT4 in T cells could contribute to allograft failure. By applying LASSO logistic regression on the randomly selected samples from GSE21374, which were regarded as the training cohorts, seven potential targets (CD247, NKG7, CD6, CCL5, FGD3, APBBAIP, and ARHGAP30) from the hub gene sets were determined and used to establish a diagnostic model ( Supplementary Figure 7A; Supplementary Table 5 ). The diagnostic ability of the model was further tested in the training cohorts ( Supplementary Figures 7A–C ), and in the rest of the samples of GSE21374, which was determined as the validation cohort ( Figure 3G ), of which the ROC curves and the overall survival analysis showed a high diagnostic and prognostic power of the model. Among the seven genes, CD247, NKG7, CD6, and CCL5 were presented relevant to transplantation rejection, reflecting poorer allograft survival after renal transplantation (49–52). Intriguingly, all seven genes were remarkably expressed in HR especially CD6, CD247, NKG7, and CCL5, which were specifically expressed in T cells ( Figures 3H, I ). These results revealed that STAT4 as a core transcription factor, mediated T-cell immune responses, which is essential in HR progression and renal allograft failure.

Since STAT4 and almost all of its putatively target genes triggered adverse allograft survival in patients who received renal transplantation, the essential upstream regulators of STAT4 deserve further identification. To explore upstream signaling pathways targeting STAT4, Nichenet’s (36) analysis was first performed to determine the overall ligand and receptor pairs targeting STAT4 in T cells from the HR group ( Figure 4A ). All of the receptors, signaling mediators, and transcription factors (TFs) in the network were picked and those uncorrelated to HR were filtered out by limiting the p-value larger than 0.05, of which IL6ST, MET, and CXCR4 were verified as upstream signaling molecules to regulate STAT4 (53–55) ( Figure 4B ). We applied these signaling molecules to the LASSO logistic regression on the training cohorts, which were randomly selected from GSE21374, and eight optimal genes, namely, CD44, FTH1, CXCR4, PTPN6, PRDX2, EWSR1, UBB, and RPS19BP1, were employed to establish a diagnostic model ( Figure 4C ; Supplementary Table 5 ). ROC curves and overall survival analysis revealed a high diagnostic and prognostic power of the model on both the training cohorts ( Supplementary Figure 8A ) and the validation cohorts, which consisted of samples in GSE21374 excluding those in the training set ( Figure 4D ). Among these genes, CD44, CXCR4, PRDX2, and UBB were significantly related to transplantation rejection and poor survival ( Figure 4E ), which were also proved by researchers’ studies (56–59). The other four were newly discovered genes potentially playing key roles in rejection and graft failure ( Figure 4E ). Interestingly, we found that CD44, CXCR4, PTPN6, and EWSR1 were specifically expressed in an HR state from the new classifications in bulk RNA-seq ( Figure 4F ). As expected, the overall survival probability of patients with higher expression of CD44, CXCR4, and PTPN6 showed worse disease consequences ( Figure 4G ), suggesting that these genes were HR regulators of allograft loss. It is known that antagonists of CD44 and CXCR4 can help improve outcomes in allograft rejection (56, 60), which further supports our hypothesis.

Importantly, PTPN6, which was newly identified by our gene regulatory network, encodes a member of the protein tyrosine phosphatase (PTP) family and regulates multiple cellular processes, such as cell growth, differentiation, and tumorigenesis (61). It is demonstrated that PTPN6 functions in a TCR-dependent manner (62) and elevated expression of PTPN6 recruits infiltration of T cells (63), illustrating its role in mediating STAT4-induced T-cell immune responses.

To confirm the precision of the unsupervised classification built in our study, we applied this strategy to the external microarray dataset (GSE21374) and yielded 4 rejection states ( Figure 5A ). Each status was finely characterized and annotated by the signature gene sets, including STA (v-STA, NECTIN1 ⁺), Fib (v-Fib, CA3 ⁺), Prog2 (v-Prog2, SLC5A3 ⁺), and HR of validation (v-HR, STAT4 ⁺). Therefore, both the method and signature genes identified in rejection states worked well in these external data. Importantly, PTPN6, as well as CD44 and CXCR4, were upregulated specifically in HR in the validation cohort, demonstrating that the two known regulators, especially the newly discovered gene PTPN6, were critical in HR progression ( Figure 5B ).

To further confirm the HR characteristic prone to allograft loss in v-HR, we isolated samples with a final diagnosis of allograft loss, and divided them into non-HR (not belonging to v-HR) and v-HR according to the new classification. Indeed, survival probability analysis revealed that patients from v-HR showed greater susceptibility to renal failure compared to the non-HR group ( Figure 5C ).

Taken together, our new classification strategy proves to be precisely helpful for distinguishing renal allograft rejection status, upon which we also provide new insights into PTPN6-STAT4-immune responsive signaling in T cells that mediates HR rejection ( Figure 5D ).