Spatial atlas of diabetic kidney disease reveals a B cell-rich subgroup

Sample acquisition and ethics approval

The University of Pennsylvania institutional review board (IRB) approved the collection of human kidney tissue for this study. Left-over kidney samples were irreversibly de-identified, and no personal identifiers were gathered. Therefore, they were exempt from IRB review (category 4). We engaged an external, honest broker who was responsible for clinical data collection without disclosing personally identifiable information. Participants were not compensated. Some samples were collected as part of the Transformative Research in Diabetic Nephropathy (TRIDENT) study, a multicentre observational cohort study that enrolls, follows, and performs multiomics characterization of individuals with DKD. This study was approved by the University of Pennsylvania IRB and by the TRIDENT Steering Committee. Informed consent was obtained from each participant. The present research using the UKBB Resource was approved under Application Number 273810. Participants from the UKBB provided written informed consent allowing the use of their samples and data for medical research purposes.

CosMx sample preparation and data preprocessing

Tissue sections were cut at 5 µm thickness and prepared according to the manufacturer specifications (NanoString Technologies). We used the human universal cell characterization RNA probes, and a 50 additional custom probes for the following genes: ESRRB, SLC12A1, UMOD, CD247, SLC8A1, SNTG1, SLC12A3, TRPM6, ACSL4, SCN2A, SATB2, STOX2, EMCN, MEIS2, SEMA3A, PLVAP, NEGR1, SERPINE1, CSMD1, SLC26A7, SLC22A7, SLC4A9, SLC26A4, CREB5, HAVCR1, REN, AP1S3, LAMA3, NOS1, PAPPA2, SYNPO2, RET, LHX1, SIX2, CITED1, WNT9B, AQP2, SCNN1G, ALDH1A2, CFH, NTRK3, WT1, NPHS2, PTPRQ, CUBN, LRP2, SLC13A3, ACSM2B, SLC4A4, PARD3, XIST and UTY. We used the DAPI, CD298/B2M, CK8/18 and PanCK/CD45 for additional staining per Nanostring protocol. Imaging was performed using configuration A. After imaging was completed, the flowcell was incubated in 100% xylene overnight, the coverslip was removed from the slide with a razor blade—the slide was then stained with haematoxylin and eosin. The expression matrix and metadata from each CosMx run was exported from the AtoMx platform and was converted to a Python object using squidpy. All samples were merged and preprocessed and analysed together using scanpy. Cells with fewer than 30 counts were filtered out.

Xenium sample preparation and data preprocessing

Tissue sections were cut at 5 µm thickness and cut onto a Xenium slide. according to the manufacturer specifications (10X Genomics). We used the human Xenium Prime 5K Human Pan Tissue and Pathways Panel with 100 additional custom probes for the following genes: TPM1, ESRRB, COL6A3, AGR2, SLC26A7, ATP1B1, SLC8A1, ATP6AP2, TAGLN, SPP1, SAT1, MYL9, LDB2, DEFB1, COL1A2, ACTA2, ST6GALNAC3, SLC13A3, SLC12A3, SLC12A1, MGP, IGHG1, FN1, C7, ACSM2B, AIF1, APOE, AQP3, AZGP1, C1QA, C1QB, C1QC, CAV1, PPIA, CD74, CHI3L1, COL1A1, COL6A1, CRYAB, CXCL14, ENO1, HLA-DPA1, HLA-DRA, IFI27, IGHA1, IL32, KLF2, LGALS3, LUM, MMP7, PIGR, S100A2, SLC4A4, SLPI, SOD2, SPINK1, SOX4, SPOCK2, TACSTD2, TM4SF1, TPM2, VIM, ZFP36, AQP2, RNASE1, ALDOB, PCGF6, RHOB, CD81, ASS1, MYL6, COX8A, CTSB, GATM, MT1G, TMSB10, COL3A1, MIF, TPT1, COL6A2, BST2, CLU, APOC1, APOD, PHKG2, RGCC, HLA-DQA2, CORO1A, HSPB1, ADIRF, CKB, HLA-DQB1, COX5B, MT1H, RAMP3, TYROBP, LAMTOR5, ITM2B, UBB and CTSD. Additionally, the same sections were stained according to the Xenium Cell Segmentation workflow for automated morphology-based cell segmentation, and subsequently loaded onto the Xenium Analyzer for in situ transcriptomic analysis. Xenium raw output files were processed using the spatialdata framework (v0.0.14) with the spatialdata-io Xenium plugin. Xenium transcript and segmentation data were loaded from the manufacturer’s output directory using the xenium() function, which parses transcript tables, cell segmentation boundaries and spatial metadata into a structured SpatialData object. The gene expression table was extracted as an AnnData object for downstream single-cell analysis. Cells with fewer than 30 counts were filtered out.

Integration of the single-cell spatial kidney atlas

To integrate snRNA-seq with spatial transcriptomic data, we applied a deep generative modelling framework using scVI and scANVI, as implemented in the scvi-tools Python package (v1.0+). Two parallel integrations were performed: one with CosMx and one with Xenium data. Raw count matrices were used, and cells were annotated with metadata, including technology, project, and sample identifiers. The SCVI.setup_anndata() function was used to register these batch and covariate annotations. The scVI model was initialized with a negative binomial likelihood and trained using 4 hidden layers and a 30-dimensional latent space. We applied early stopping based on validation loss and used a batch size of 2,048. We refined this by training a scANVI model based using the scVI-trained model as a prior and using the annotations from the annotated snRNA-seq reference atlas. scANVI was trained with early stopping, The mean distance between a given CosMx cells and 10 snRNA-seq cells was then calculated in the resulting latent space using k-nearest neighbours. Cells with a high distance were discarded and reintegration was performed. To create a unified reference map across technologies, we then concatenated the filtered CosMx and Xenium objects. scVI was trained using the raw counts as input, with tech as batch key, and proj, orig_ident, nCount_RNA, and percent_mt as additional covariates. The model was trained using 30 latent dimensions, 4 hidden layers, and a negative binomial likelihood, with dropout and early stopping enabled. scANVI was subsequently initialized using the trained scVI model and fine-tuned using coarse annotations. The model was trained for 50 epochs with a learning rate of 0.001. UMAP coordinates were computed based on the scANVI latent space, and Leiden clustering was performed for annotation. Based on the final scanVI latent space, the distance from each spatial cell was calculated to the nearest snRNA-seq cell using a k-nearest neighbour graph based on Euclidean distance in the latent space. The expression of 3,000 highly variable genes of the snRNA-seq dataset was then inferred for each spatial cell by averaging the corresponding gene expression values of its nearest neighbours, with contributions weighted by the squared inverse of their latent space distances, such that closer neighbours contributed more strongly.

Annotation-based niche generation

To quantify cell-type-specific neighbourhood compositions across multiple spatial scales, we computed the number of neighbouring cells of each cell type within 20, 40, 60 and 80 μm radii for each cell in the combined CosMx and Xenium dataset. The input contained spatial connectivity matrices at each distance threshold. For each sample, we extracted subsets of the data and iteratively computed, for every cell type, the total number of neighbours per cell using the respective spatial connectivity matrix. This process yielded four neighbourhood matrices, each representing a different spatial scale. To capture neighbours newly included at increasing distances, we computed differences between consecutive shells (for example, 40 μm minus 20 μm). Cells with zero neighbours across all types were excluded, and the bottom 5% of sparse neighbourhoods (based on total neighbour count at 80 µm) were filtered out to mitigate edge effect. Cell-type-specific neighbour compositions across all four scales were then used as input features to define spatially resolved kidney niches. These features were normalized using StandardScaler, and dimensionality reduction was performed using principal component analysis (PCA). The optimal number of clusters was determined using the elbow method, followed by MiniBatchKMeans clustering to assign niche labels. To assess disease-associated shifts in kidney niche composition, we quantified the relative abundance of annotated niche labels across healthy and DKD samples. For each condition and niche, we calculated the number of contributing cells and normalized this count by the total number of cells from the respective condition to correct for differences in sampling depth. This yielded a normalized abundance matrix for each niche within control and DKD groups. To compare niche representation between groups, we computed the relative increase or depletion in DKD by calculating the ratio of normalized counts (DKD/control), assigning negative values when niches were more abundant in controls. These relative changes were visualized as a heat map.

COVET and Hotspot analysis

We used the covariance environment (COVET) methodology²¹ on all cells in the CosMx dataset with more than 50 neighbours within a 200 µm radius. This approach incorporates the gene–gene covariate structure with spatial location considerations. Inputs to the COVET matrix included: (1) spatially resolved gene expression; (2) 64 highly variable genes (selected via Scanpy’s highly_variable_genes); and (3) the number of spatial neighbours (k = 8). For each cell in the CosMx dataset with at least 50 neighbours within 200 μm, we computed an adjusted covariance matrix across the 64 genes to capture niche-specific co-expression patterns. Pairwise niche distances were then derived using approximate optimal transport (AOT). Based on this distance matrix, we calculated UMAP for visualization after PCA which help to better examine the relationships within the COVET niche representation. For analysis of spatially correlated transcripts, we used Hotspot²³ for each niche, based on the niche-specific latent representations calculated by COVET analysis. Specifically, we first construct a cell–cell similarity graph in which nodes represent cells and edges encode similarity in the COVET latent space. Gene–gene colocalization is then quantified using a local correlation statistic that evaluates whether two genes tend to be expressed in the same regions of this similarity graph. For a pair of genes \(x\) and \(y\), the local correlation is defined as:

where \(\tildex_i\) and \(\tildey_i\) denote standardized expression values of genes \(x\) and \(y\) in cell \(i\), respectively, and the sum is taken over all neighbouring cell pairs \((i|j)\) in the similarity graph with edge weights \(w_ij\). The difference in gene–gene correlations between DKD and control was calculated by subtracting the control matrix from the DKD matrix. To normalize values across datasets, each correlation matrix (control, DKD and the difference matrix) was scaled to the range [–1, 1] using min–max scaling. Hierarchical clustering was applied to the scaled control matrix to define a consistent gene ordering, which was subsequently used to reorder all matrices for visualization. Clustered heat maps were generated using Seaborn’s clustermap and heat map functions, with selected gene pairs highlighted to visualize conserved and differential co-regulation.

We visualized gene pairs with absolute correlation values above 0.99 in both control and DKD, identifying robust co-expressed gene modules. Additionally, we computed the top 30 gene pairs with the greatest correlation differences between DKD and control and visualized these.

NicheCompass analysis

To infer spatially defined gene programmes and cellular niches, we trained the NicheCompass model on all cells in the CosMx dataset with more than 50 neighbours within a 200 µm radius. The model was configured with raw count data, sample-specific categorical covariates, and a spatial adjacency matrix for edge-aware graph neural network learning. The encoder used a graph attention convolutional layer (gatv2conv), and categorical covariates were injected into the decoder layer. Gene programme activity thresholds were determined using a relative cutoff (active_gp_thresh_ratio=0.01). The model was trained for 50 epochs and 25 additional epochs across all gene programmes, using the following optimization parameters: learning rate (lr = 0.001), edge reconstruction penalty (λ_edge_recon = 5 × 10⁵), gene expression reconstruction penalty (λ_gene_expr_recon = 300), and L1 regularization for novel gene programmes (λ_L1_addon = 30). The batch size for edge sampling was set to 1,900. We computed a neighbourhood graph and UMAP embedding based on the latent representation and performed clustering using the Leiden algorithm. Clusters were manually annotated to define discrete cellular. To identify gene programmes enriched in specific niches, we performed differential gene programme testing across latent clusters using a log Bayes factor threshold of 2. The resulting active gene programmes were visualized via heat maps after min–max normalization, highlighting niche-specific transcriptional signatures.

Injured tubular microenvironment generation

To characterize injured tubular microenvironments, we subset the CosMx and Xenium datasets to include injured proximal and thick ascending limb cells (iPT, iTAL). For each of these cells, we quantified the surrounding neighbourhood composition across four spatial shells (20 μm, 40 μm, 60 μm and 80 μm) and aggregated related cell types (for example, glomerular, collecting duct and endothelial) to reduce dimensionality. The resulting neighbourhood matrices were normalized using StandardScaler, and dimensionality reduction was performed using PCA. We applied MiniBatchKMeans clustering to identify niche subtypes. Clusters were manually annotated based on dominant neighbourhood composition and labelled as profibrotic, perivascular, periglomerular, proximal or distal microenvironments. We then examined the cell type composition of each niche and visualized sample-specific niche distributions. Differential gene expression analysis was performed within each injured tubular subtype (iPT and iTAL), stratified by microenvironment, across both CosMx and Xenium platforms. To quantify disease-associated shifts in injured tubular microenvironments, we compared the frequency of iPT/iTAL microenvironment labels between DKD and control samples. For each sample, the number of cells assigned to each microenvironment was calculated and normalized to the total number of PT, iPT, TAL, or iTAL to correct for sampling differences. We compared normalized microenvironment frequencies between DKD and control samples, computing relative changes as fold increases (DKD/control) or depletions (−control/DKD). The resulting values were visualized in a heat map.

Immune microenvironment generation

To construct spatial neighbourhood profiles centred around immune cells, we used the immune cell type annotation from the kidney immune cell atlas. For each sample in the integrated Xenium and CosMx dataset, we computed cell-type-resolved neighbourhood matrices at four spatial scales (20 µm, 40 µm, 60 µm, 80 µm) using precomputed spatial connectivity graphs. For each radius, we quantified the number of neighbouring cells belonging to each cell type per index cell, resulting in four DataFrames capturing neighbourhood composition. To isolate newly added neighbours at increasing spatial distances, we calculated the difference between consecutive shells (e.g., 40 µm minus 20 µm). Cells with zero neighbours across all types were excluded, and the bottom 5% of sparse neighbourhoods (based on total neighbour count at 80 µm) were filtered out to ensure robust representation. Each DataFrame was then normalized by the total number of neighbours per cell at that distance to yield cell-type-specific neighbour fractions, and column names were suffixed with the corresponding spatial scale. Finally, cells shared across all four scales were retained and the resulting matrices were concatenated into a single comprehensive neighbourhood profile matrix. This matrix served as the input for downstream immune microenvironment classification. Using these precomputed neighbour fraction matrices at 20–80 µm radii, we aggregated cell types into broader tissue compartments, such as glomerular cells, injured tubules, collecting ducts, loop of Henle, endothelium, and healthy tubules. The resulting DataFrame was standardized and subjected to PCA. MiniBatchKMeans clustering was then performed on the PCA-reduced data. The optimal number of clusters was determined using the elbow method. Cells were assigned cluster labels corresponding to immune microenvironment types. To visualize immune cell composition within each immune microenvironment, we grouped cells by immune microenvironment and immune cell type, and counted the number of occurrences to obtain cell type frequencies per microenvironment. To simplify visualization, rarer immune cell types (basophils/mast cells, conventional and plasmacytoid dendritic cells, NK cells, neutrophils) were aggregated into a single “other” category. We then generated pie charts representing the immune composition of each microenvironment using Matplotlib. To visualize the full cellular composition of the 20 µm immune microenvironments, we grouped cells in the 20 µm surrounding each immune microenvironment and cell type, then counted occurrences to calculate the local composition. To reduce complexity and highlight key structures, related cell types were collapsed into broader categories: PT and TAL into “PT_TAL”, iPT and iTAL into “iPT_iTAL”, glomerular cells (Podo, MC1, PEC, EC_glom) into “glom”, and vascular populations (EC_Peritub, EC_DVR, VSMC) into “ec”. Remaining cell types (e.g., DCT, CNT, intercalated cells) were grouped into an “other” category. We then generated pie charts representing the cellular composition of each microenvironment using Matplotlib.

To assess changes in immune cell composition within the glomerular microenvironment, we calculated the frequency of immune cell subtypes by counting the number of cells per immune subtype and normalizing to the total number of cells in the glomerular immune microenvironment in control and DKD samples. Samples with fewer than 50 glomerular cells were excluded. Fractions were centred by subtracting the median of control samples for each immune subtype. Group-wise comparisons between DKD and control were performed using Welch’s t-test. Results were visualized using box plots of centred immune cell fractions. To analyse immune cell composition in the injured tubular immune microenvironment, we calculated the frequency of immune cell subtypes by counting the number of cells per immune subtype and normalizing to the total number of cells in the injured tubular immune microenvironment in control and DKD samples. Immune cell subtype counts were computed per sample and normalized by the total number of cells in the corresponding niche to obtain cell-type fractions. Samples with fewer than 50 total cells in the injured tubular immune niche were excluded. Immune cell fractions were centred by subtracting the median fraction of control samples. Statistical comparisons between DKD and control samples were performed using Welch’s t-test. Centred fractions were visualized using violin and box plots to assess relative enrichment or depletion of specific immune cell types within the tubular immune microenvironment in DKD.

Human kidney immune cell atlas generation

For the integration of the human kidney immune cell atlas, we integrated immune cells from spatial transcriptomic data with immune cells from snRNA-seq, we used a deep generative modelling framework based on scVI and scANVI, implemented in the scvi-tools Python package (v1.0+). Immune cells were first subset separately from the CosMx and Xenium spatial datasets. Each subset was integrated independently with immune cells from the snRNA-seq reference atlas. Raw count matrices were used and additional covariates included technology, project, and sample identifiers were used. The SCVI.setup_anndata() function was used to register batch and covariate annotations. scVI was trained with a negative binomial likelihood, four hidden layers, and a 30-dimensional latent space and early stopping. After training, scANVI was initialized using the scVI model and refined with annotations from the snRNA-seq reference. scANVI was trained with a learning rate of 0.001. We then trained a scANVI model using known annotations and predicted cell identities for Xenium/CosMx cells labelled as ‘Unknown’. To refine predictions, we iteratively reassigned Xenium/CosMx cell types based on their 10 nearest neighbours in the scANVI latent space, using majority voting if their mean distance fell below a defined percentile threshold (initially 50%, increased to 70% and 90% in subsequent iterations). After each round, updated labels were used to retrain the scANVI model, leading to progressively improved annotation accuracy. Following independent integration of CosMx and Xenium immune cells, the resulting AnnData objects were concatenated. A final scVI model was then trained on the combined dataset using raw counts, with tech as the batch key and proj, orig_ident, nCount_RNA, and percent_mt as covariates. scANVI was reinitialized on the combined latent space and refined using coarse annotations. UMAP coordinates were computed from the final scANVI latent space, and Leiden clustering was applied. Expression of immune cell genes was then inferred for each spatial cell by averaging the corresponding gene expression values of its nearest neighbours, with contributions weighted by the squared inverse of their latent space distances, such that closer neighbours contributed more strongly.

B cell atlas generation

Raw count matrices for B cells and plasma cells were provided in Matrix Market Exchange (MTX) format from a recently published single-cell B cell atlas³⁸, along with matching cell-level metadata (CSV) and gene name annotations. The count matrices were loaded using scipy.io.mmread and transposed to ensure cells corresponded to rows and genes to columns. Each count matrix was then combined with its corresponding metadata and gene names to construct a sparse AnnData object using Scanpy (v1.9+). To construct a unified B cell atlas, we first excluded a curated list of blacklisted genes prior to downstream modelling. The dataset was normalized using total count normalization followed by log-transformation. We identified 3,000 highly variable genes using the Seurat flavour in Scanpy with dataid as the batch key. Cells with fewer than 200 counts were removed. The resulting dataset was used as input for the scVI model. We initialized a deep generative model using scVI (v1.0+, scvi-tools) with the raw count layer as input, dataid as batch, and additional covariates including orig_ident, nCount_RNA, and percent_mt. The model used a negative binomial likelihood, three hidden layers, and a 30-dimensional latent space. Training was performed with a learning rate of 0.001, early stopping. For cell type refinement, we used scANVI and coarse labels from the original annotation. The scANVI model was initialized from the trained scVI model with celltype_coarse as the label key and Unknown as the unlabelled category. It was trained for 100 epochs using a learning rate of 0.001, batch size of 256, and n_samples_per_label = 50, with early stopping enabled. Leiden clustering was applied to the scANVI latent space and cell types assigned to known B cell states including BIe, atypical memory B, age-associated B cells, germinal centre B cells and plasma cells. To enable integrated analysis of single-cell and spatial transcriptomic B cell data, we concatenated annotated single-cell RNA-seq data with B and plasma cells extracted from Xenium datasets. To transfer annotations from single-cell RNA-seq to Xenium spatial transcriptomics data, we trained using scVI and scANVI. Annotated single-cell and unannotated Xenium B and plasma cells were concatenated and filtered to include only samples with at least 50 cells. The scVI model was trained using tech as the batch key and proj and orig_ident as categorical covariates, with nCount_RNA and percent_mt as continuous covariates. We then trained a scANVI model using known annotations and predicted cell identities for Xenium cells labelled as Unknown. To refine predictions, we iteratively reassigned Xenium cell types based on their 10 nearest neighbours in the scANVI latent space, using majority voting if their mean distance fell below a defined percentile threshold (initially 50%, increased to 70% and 90% in subsequent iterations). After each round, updated labels were used to retrain the scANVI model, leading to progressively improved annotation accuracy. Following scVI/scANVI integration, we computed UMAP embeddings based on the scANVI latent space, and Leiden clustering was applied. Based on the clustering results and expression profiles, cluster identities were manually curated and mapped to B cell subtypes, including naive B cells, memory B cells, activated B cells, germinal centre B cells, atypical memory and plasma cells.

Imaging mass spectrometry

FFPE tissue sections were first deparaffinized and rehydrated following standard protocols. Antigen retrieval was carried out using citrate buffer (pH 6.0), followed by blocking with 5% BSA. Primary antibodies targeting proteins of interest were conjugated to cadmium (Cd) metal isotopes using the Maxpar Antibody Labeling Kit (Fluidigm) and incubated with the tissue sections overnight at 4 °C. After washing, metal-conjugated secondary antibodies were applied to enhance signal amplification. Nuclear staining was performed using an iridium intercalator. IMC was conducted on a Hyperion Imaging System (Fluidigm) equipped with a 20× objective. Tissue ablation was performed using a laser with 1 μm resolution, and metal isotope signals were quantified using a time-of-flight mass spectrometer.

Differentially expressed genes

Raw counts were normalized, and log-transformed DEGs were calculated using Wilcoxon rank-sum test with the built in scanpy function. DEGs between patient subcluster in the bulk RNA-seq of the TRIDENT cohort were calculated using Limma⁶⁵.

Pathway enrichment analysis using Enrichr

To identify pathways enriched in niche-specific upregulated genes, we performed enrichment analysis using the gseapy Python package (v1.0.6)⁶⁶. Differential gene expression results for each niche were loaded from precomputed.csv files, and low-confidence genes with adjusted P values of zero were replaced with a small non-zero constant (10⁻³¹⁰). Only genes with positive log fold change and an adjusted P value < 0.05 were retained for enrichment. The background gene set was defined as the union of all tested genes across niches. Gene set enrichment was performed using the enrichr function with the WikiPathways 2023 Human gene set collection (min size = 20). For each niche, significantly enriched pathways (P < 0.05) were extracted from the Enrichr results. Enrichment results were visualized using dot plots or stored for downstream comparison.

Pseudobulk differential expression analysis and spatial gene signature derivation

To identify niche-specific gene signatures, we performed pseudobulk differential expression analysis using the CosMx dataset. We aggregated raw gene expression counts by summing expression values across all cells from the same sample and annotated niche, resulting in a pseudobulk expression matrix (cells grouped by sample and niche). Corresponding metadata containing sample and niche labels was generated. The resulting count matrix and metadata were exported and imported into R for differential gene expression analysis. We then used the DESeq2 package to perform differential expression for each annotated niche versus all other niches. For each contrast, we created a binary group label (‘niche’ versus ‘non-niche’) and ran the DESeq2 pipeline with estimateSizeFactors (type = “poscounts”) and DESeq. Results were extracted using the results() function. DEGs were defined as those with an adjusted P value < 0.05 and absolute log₂ fold change >0.5. To quantify significance and ranking, a π-score was computed as the signed −log₁₀ of the adjusted P value. Genes were filtered for adjusted P value < 0.05 and log₂ fold change >0.5. To rank genes by both significance and magnitude of change, we computed a π-score for each gene, defined as the product of the sign of the log₂ fold change and the negative log₁₀ of the adjusted P value. This composite metric prioritizes genes with strong upregulation and high statistical confidence. The top 10 upregulated genes per niche were selected as spatial gene signatures.

Ligand–receptor interactions

To identify microenvironment-specific cell–cell communication signals, we performed the statistical ligand–receptor interaction analysis using CellPhoneDB v5 (ref. ⁶⁷. An AnnData object containing imputed gene expression was used as input, along with metadata specifying cell-type annotations and microenvironment labels. Cells were first filtered to exclude minor cell types comprising less than 1% of cells within each niche. For downstream compatibility with CellPhoneDB, we annotated each cell with its corresponding microenvironment and defined composite labels combining microenvironment and cell type. These metadata were formatted into two CSV files: one mapping barcodes to cell types and another associating each cell with a microenvironment. We used the CellPhoneDB statistical framework with 1,000 iterations and a threshold of 1% gene expression to determine significant ligand–receptor pairs. Interactions were scored, and significance was determined using a P value cutoff of 0.05. The HGNC gene symbol was used for annotation, and the analysis was performed using 8 computational threads. Results were saved for downstream interpretation and visualization. For calculation of ligand–receptor interactions, we used the imputed gene expression. To identify disease-relevant cell–cell interactions associated with kidney function, we used CellPhoneDB statistical output tables. Interaction scores were indexed by ligand–receptor pairs and filtered to include only interactions between specific sender and receiver cell types within the glomerular niche. For each sender–receiver pair (for example, EC_glom to Podo), we calculated Pearson correlations between interaction scores and GFR across control and DKD samples. Only interactions expressed in more than eight samples were retained for robust correlation estimates.

Bulk RNA-seq correlation

Processing was performed as previously described²⁸. In summary, quality control of the sequencing data was assessed using FASTQC. Adaptors and low-quality bases were removed with TrimGalore (v0.4.5). The trimmed FASTQ files were then aligned to the human genome (hg19/GRCh37) using STAR (v2.7.3a) based on Genecode v19 annotations. Gene expression levels were quantified with RSEM by calculating uniquely mapped reads as transcripts per million (TPM). Count matrices were log-transformed, and the indicated genes extracted. Gene expression values were centred and min–max normalized (range: –1 to 1). For each gene, Pearson correlation was computed. Scatterplots with regression lines and correlation coefficients were generated using Seaborn. For each spatial gene signature, the corresponding 10 genes were extracted. Each gene was scaled to a range from −1 to 1, and the scaled values of the 10 genes within each spatial gene signature were summed. The resulting score was then correlated with the fibrosis score of the respective sample using Pearson correlation.

Hierarchical clustering of bulk RNA kidney samples

Hierarchical clustering was conducted on the scaled TPM matrices using 10 genes from the B predominant immune microenvironment gene list as a reference (Supplementary Table 23). The datasets were clustered using Ward’s method with Euclidean distances.

TRIDENT plasma proteomics analysis

We trained an elastic-net logistic regression classifier using tidymodels⁶⁸. Because B cell classes were imbalanced, we applied the synthetic minority oversampling technique (SMOTE) to the training data. Model complexity (α, λ) was tuned via fivefold cross-validation to maximize the area under the receiver operating characteristic curve (AUROC). Predictors with non-zero coefficients in the refit model were then used as covariates in a Cox proportional hazards model of kidney survival, defined as a composite outcome of death, dialysis, or ≥40% decline in eGFR. We compared the C-statistic of a model including clinical characteristics plus these proteins with a clinical-only model. To account for potential overfitting, we performed 1,000 bootstrap resamples and reported the mean C-statistic for each model.

Validation in UK Biobank cohort

The UK Biobank (UKBB) is a large-scale, community-based cohort comprising over 500,000 participants aged 40–69 years at recruitment (2006–2010) from 22 assessment centres across the UK. For validation purposes, we included 3,309 participants who underwent plasma proteomic profiling and were diagnosed with diabetes mellitus¹. This subset is representative of the overall UKBB population. Primary outcomes in the UKBB included end-stage kidney disease (ESKD) (identified via ICD-10 codes N185 and N180), dialysis or kidney transplantation. ESKD diagnoses and dates were determined using first-occurrence records from multiple data sources. eGFR was calculated using the CKD-EPI formula based on serum creatinine and cystatin C measurements obtained from baseline and follow-up records. (Study Application 273810). The predictive performance of the developed models was evaluated using time-dependent receiver operating characteristic curve (tAUC) were calculated to assess the model’s performance at different time points⁶⁹.

Immune complex pull-down, LC–MS acquisition and data analysis

Kidney tissues from nine controls and nine samples with B cell-predominant inflammation from the Susztaklab Biobank were lysed under native conditions to preserve protein–protein interactions. Immunoglobulin and associated proteins were captured using a biotin-conjugated anti-IgA antibody (ThermoFisher, 7102882100) and pulled down with Dynabeads MyOne Streptavidin C1 (ThermoFisher, 65001). The bead-bound complexes were prepared for liquid chromatography–mass spectrometry (LC-MS) analysis following an on-bead digestion workflow: samples were reduced with tris(2-carboxyethyl)phosphine (TCEP), alkylated with iodoacetamide (IAA) and digested overnight with a Trypsin/LysC mix (Promega, V5072). The resulting peptides were then labelled with tandem mass tags (TMTpro; ThermoFisher, A44520 and A52045), multiplexed, and cleaned prior to high pH reversed-phase peptide fractionation to achieve deep proteome coverage. The resulting fractions were analysed on an Orbitrap Eclipse mass spectrometer equipped with FAIMS Duo and real-time search for in-depth quantitative profiling. Raw data were processed using the GFY core pipeline and statistically analysed with the msTrawler framework⁷⁰. Differential expression and covariation analyses were performed to identify patterns of protein interactions, and IgA correlation metrics were combined to generate a ranked list of IgA-interacting proteins.

Statistics and reproducibility

Statistical analyses are described in the Methods and figure legends. All tests were two-sided unless stated otherwise. Differential expression was assessed using Wilcoxon or limma with Benjamini–Hochberg correction. Group comparisons used Welch’s t-test, survival differences were evaluated by log-rank test, and associations were assessed by Pearson correlation with 95% confidence intervals. Pathway enrichment was performed using GSEApy (GSEA or over-representation analysis with Fisher’s exact test and Benjamini–Hochberg correction). No statistical method was used to predetermine sample size. The spatial transcriptomics discovery cohort comprised 64 human kidney samples (48 CosMx, 16 Xenium), and the snRNA-seq cohort included 150 samples. Sample size was determined by tissue availability and feasibility of high-resolution spatial profiling. The study was observational and not randomized. Samples were assigned to groups based on established clinical diagnoses. Investigators were not blinded during analysis; however, spatial transcriptomics data processing were conducted using automated computational pipelines, minimizing subjective bias. Representative histological and spatial images reflect patterns observed in the indicated samples. Biological replicates correspond to patient samples; in cases where samples were profiled twice, measurements were treated as independent technical datapoints. Sample-level comparisons were performed on aggregated per sample values unless otherwise stated. Samples and cells were excluded based on predefined quality control criteria (see Methods). Robustness was supported by cross-platform replication (CosMx and Xenium), integration with snRNA-seq data, external validation cohorts and cross-validation or bootstrap resampling where applicable.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.