Pathology-oriented multiplexing enables integrative disease mapping

Archival tissues

Human samples

FFPE tissues were collected and prepared according to institutional protocols. For Fig. 2, validation kidney biopsy specimens from individuals with ANCA-associated CGN were obtained from the Hamburg Glomerulonephritis Registry (https://www.sfb1192.de/en/register). For Fig. 3, control kidney specimens were obtained from nephrectomies performed on individuals with renal cell carcinoma in collaboration with the Division of Nephrology and Clinical Immunology, RWTH Aachen University Medical Center. Kidney biopsy samples from individuals with DKD were obtained from the Department of Nephrology and the Department of Pathology Georges Pompidou European Hospital, Assistance Publique–Hôpitaux de Paris. For Fig. 4, nested protocol research kidney biopsy samples were obtained from volunteers (adolescents and young adults; n = 13) with T2D (12–21 years of age, T2D onset at <18 years of age, T2D duration 1–10 years and HbA1c < 11%) from the Renal HEIR and the IMPROVE-T2D studies. The participants were recruited from the Type 2 Diabetes and Metabolic Bariatric Surgery clinics at the Children’s Hospital Colorado Anschutz Medical Campus in Aurora. T2D was defined according to criteria of the American Diabetes Association plus the absence of glutamic acid decarboxylase, islet cell, zinc transporter 8 and/or insulin autoantibodies. The Renal HEIR and IMPROVE-T2D cohorts have intentionally harmonized study protocols. Medication use was recorded for all participants, and T2D treatment, including SGLT2 inhibitors, was determined by their medical provider. Normative kidney reference tissue from research biopsy samples were provided by five healthy young adult participants in the CROCODILE study (NCT04074668). For Supplementary Figs. 3 and 4, kidney biopsy samples were obtained from the Hamburg Glomerulonephritis Registry (https://www.sfb1192.de/en/register), liver specimens were provided by the Institute of Pathology, University Medical Center Hamburg-Eppendorf and brain specimens were provided by the Institute of Neuropathology, Freiburg University Hospital. Ethics approvals were obtained from the Institutional Review Board of the RWTH Aachen University Medical Center (EK-016/17), the local ethics committees of the Chamber of Physicians in Hamburg (PV4806) and Freiburg (Ethikvotum 10008/09), the Paris Ethics Committee (IRB00003888, FWA00005831) and the Colorado Ethics Committee (NCT03584217 and NCT03620773). All tissue collections were performed in accordance with the ethical principles stated by the Declaration of Helsinki.

Rodent samples

Archival FFPE tissues from experimental immune-mediated kidney disease and DKD were collected according to institutional protocols of Hamburg, Melbourne, Heidelberg and Paris (N047/20, MMCB/2006/29, H2052-2071/23 and 358-86/609EEC, respectively). All experimental animals were housed at an ambient temperature of 20 ± 2 °C, humidity of 55 ± 10% and a light–dark cycle of 12–12 h. In brief, mouse crescentic nephritis was induced according to an established protocol⁵⁸. Rat tissues were obtained from two experimental set-ups^32,59. Administration of a JNK inhibitor (CC930, dose of 60 mg kg^–1 in 0.5% carboxymethyl cellulose) or vehicle alone was performed twice daily by oral gavage. The prevention study (therapy started at day 0 and animals were killed on day 1) was performed in outbred male Sprague–Dawley rats, as this strain is known to develop heavy proteinuria⁵⁹. The therapeutic study (therapy started at day 7 after disease induction and continued until animals were killed on day 28) was performed in inbred male Wistar Kyoto rats, which are prone to developing crescent formation. In both studies, proteinuria measurements and histopathology were performed according to standardized protocols^32,59. Btbr-Lep^ob/ob (Btbr^ob/ob) mice were obtained by crossing two heterozygous Btbr^ob/WT mice purchased from The Jackson Laboratory. This model shows morphological and physiological traits of DKD (that is, hyperglycaemia, albuminuria and glomerular hypertrophy). Wild-type littermates were used as controls.

Highly multiplexed imaging

Sample preparation

Depending on the number of samples, a suitable-sized glass surface was selected and coated with poly-d-lysine (1 mg ml^–1; Merck, A-003-E) for 30 min or with APTES (Merck, 440140) 10% v/v in acetone (Merck, 320110) for 2 min and then dried overnight before mounting the sections. We initially used poly-d-lysine for all our experiments but realized that significant lifting was progressively observed in all organs tested, including kidney, lung, colon, liver and brain. Lifting was initially mild in kidney samples but highly prominent in lung and colon specimens. For example, we observed partial but meaningful tissue lifting in 76 out of 498 ROIs (15%) by the end of 49 imaging cycles (Fig. 3). Furthermore, from 23 lung specimens analysed over 8 imaging cycles, lifting was already observed in 7 of them (30%). These observations across multiple tissue types led us to conclude that poly-d-lysine coating exhibits organ-dependent and time-dependent reliability limitations, for which we recommend potential users to perform pilot studies in their organ of interest. However, after a comprehensive literature review, we identified APTES as an ideal coating agent. Using APTES, tissue lifting occurred in only 1% of kidney specimens (Fig. 4). After coating, FFPE tissues were cut at a thickness of 2–3 μm and carefully mounted on the coated glass surface (for example, µ-Slide 2-well glass-bottom (Ibidi, 80287), µ-Slide 8-well glass-bottom (Ibidi, 80827), Cell Imaging Plate 24-well glass-bottom (Eppendorf, 0030741021) or Nexterion glass (Schott, 1868767)). To prevent dissolution of the plastic components in the chambered coverslips and plates by the solvent used for deparaffinization, the walls of each well were protected by a seal of transparent silicone (Pattex) or a ring of solvent-resistant plastic, respectively.

The following steps were performed only once before initiating the sequence of cycles.

Deparaffinization and rehydration

Samples were treated with the following set of solutions: Histo-Clear (National Diagnostic, HS-200) three cycles of 10 min each, followed by an ethanol series consisting of three cycles of 100% ethanol (10 min), two cycles of 70% ethanol (5 min), one cycle of 50% ethanol (5 min) and finally, triple immersion in double-deionized water (ddH₂O) for 5 min each.

Antigen retrieval

Samples were immersed in target retrieval solution pH 9 (Agilent, S236784-2) and heated for 15 min using a steamer (Braun; FS 3000). Afterwards they were left to cool down to room temperature for 30–60 min. Sections were then incubated for 15 min in EnVision FLEX wash buffer (Agilent, K800721-2).

Blocking

To limit nonspecific antibody binding, samples were incubated in a blocking solution consisting of 5% BSA (Merck, A7906) in Dulbecco’s PBS (Thermo Fisher Scientific, 14190094) for 1 h at room temperature. Afterwards, samples were washed three times for 5 min with wash buffer.

Elution

An elution buffer was prepared according to a previously described formulation¹⁵, which consisted of 0.5 M glycine (Carl Roth, 3908.2), 3 M urea (Merck, U5378), 3 M guanidine hydrochloride (Merck, G4505) and 70 mM TCEP (Merck, C4706) mixed in ddH₂O and adjusted to pH 2.5. Samples were incubated in elution buffer once for 5 min and then three times for 10 min on a shaker, followed by three washes of 5 min with wash buffer.

NHS-E labelling

Whenever used as a reference for alignment, NHS-E (Thermo Fisher Scientific, A10168) diluted in PBS (1:400) was added to the samples for 1 h at room temperature. After 1 h, samples were washed three times for 5 min with wash buffer.

The following steps were completed for every subsequent cycle of staining and carried out in a light-free environment to prevent the crosslinking of antibodies.

Primary antibody stain for indirect immunofluorescence

Samples were incubated with primary antibodies in EnVision FLEX antibody diluent (Agilent, K800621-2) for either 1 h at room temperature (Fig. 2) or overnight at 4 °C (Figs. 3 and 4), followed by three washes of 5 min with wash buffer. We provide confirmation of each staining pattern for every antibody in Supplementary Data 1 and 2. We also validated the practical feasibility of 1-h incubations at room temperature under non-multiplexed and multiplexed imaging conditions (Supplementary Fig. 11).

Secondary antibody and nuclear stain for indirect immunofluorescence

Appropriately matched secondary antibodies (and directly conjugated primary antibodies) and the nuclear markers DAPI (Merck, D9542, 1:200) or DRAQ5 (Abcam, ab108410, 1:200) were mixed in antibody diluent and incubated for 1 h at room temperature. Afterwards, samples were washed three times for 5 min with wash buffer.

Imaging

An imaging buffer was prepared according to a previously described formulation¹⁵, which consisted of 700 mM N-acetyl-l-cysteine (Merck; A9165) mixed in ddH₂O and adjusted to pH 7.4. Imaging buffer was added to samples for imaging and then washed three times for 5 min with wash buffer before elution.

Antibody elution

Samples were incubated with elution buffer once for 5 min and then three times for 10 min on a shaker, followed by three washes of 5 min with wash buffer.

Thereafter, these steps were repeated until the desired number of antibodies was reached. Together, each cycle (using 1 h of incubation time for primary antibodies) can be completed in under 4 h of bench work. All cycles per experiment (antibodies and order) are described in Supplementary Table 7. Periodic acid Schiff (PAS) staining was performed after the last immunofluorescence staining only in Fig. 3 following a standard protocol, including incubation with periodic acid (Th. Geyer, 3257.1) to oxidize the sections, followed by Schiff’s reagent (Merck, 1090330500) to label glycol-containing structures. The sections were then counterstained with Mayer’s haematoxylin (Agilent Technologies, S330930-2).

Primary antibodies and lectins

For human samples

ABCG2 (Santa Cruz, sc-377176; 1:200); ACE-2 (R&D Systems, AF933; 1:200); adiponectin (Thermo Fisher Scientific, MA1-054; 1:200); AIF (Cell Signaling Technology, 5318; 1:200); AKAP12 (Proteintech, 25199-1-AP; 1:600); AKR1B1 (Thermo Fisher Scientific, PA5-82915; 1:500); AKR1C1 (Thermo Fisher Scientific, MA5-32842; 1:200); alpha B-crystallin (Proteintech, 68001-1-Ig; 1:1,000); ANXA3 (Sigma-Aldrich, HPA013398; 1:200); αSMA–FITC conjugate (Sigma-Aldrich, F3777; 1:800); aquaporin 2 (Alomone Labs, AQP-002; 1:400); β-actin (Sigma-Aldrich, A5441; 1:1,500); β-catenin (Abcam, ab6302; 1:2,000); β-tubulin (Cell Signaling Technology, 2128; 1:150); calbindin-D (Sigma-Aldrich, C9848; 1:3,000); calpain small subunit 1 (Abcam, ab92333; 1:200); calpastatin (Abcam, ab244460; 1:200); calreticulin (Abcam, ab92516; 1:300); carbonic anhydrase IX (R&D Systems, AF2188; 1:50); catalase (Proteintech, 66765-1-Ig; 1:300); CD3 (Abcam, ab11089; 1:200); CD4 (R&D Systems, AF-379-NA; 1:100); CD8 (Agilent, M710301-2; 1:200); CD34 (Agilent, GA63261-2; 1:50); CD41 (Thermo Fisher Scientific, PA5-79526; 1:500); CD42b (Abcam, ab227669; 1:100); CD44 (Cell Signaling Technology, 5640S; 1:200); CD44–Alexa Fluor 647 conjugate (BioLegend, 103018; 1:200); CD68 (BioLegend, 916104; 1:200); CD79α (Agilent, M705001-2; 1:200); CD200 (R&D Systems, AF2724; 1:100); CD206 (Proteintech, 60143-1-Ig; 1:2,000); FOS (Abcam, ab190289; 1:600); claudin 1 (Abcam, ab15098; 1:500); claudin 10 (Thermo Fisher Scientific, 38-8400; 1:100); collagen I (Southern Biotech, 1310-01; 1:200); collagen III (Abcam, ab7778; 1:200); collagen IV (Abcam, ab6586; 1:200); collagen V (Abcam, ab7046; 1:100); cubilin (R&D Systems, AF3700; 1:200), cyclin B1 (Cell Signaling Technology, 12231; 1:100); cytochrome c (Abcam, ab110325; 1:200); cytokeratin 7 (Agilent, GA61961-2; 1:300); cytokeratin 8 (R&D Systems, MAB3165-SP; 1:300); cytokeratin 19 (Abcam, ab52625; 1:300); C1QA (Proteintech, 67063-1-Ig; 1:1,000); DACH1 (Sigma-Aldrich, HPA012672; 1:200); decorin (R&D Systems, AF143, 1:50); E-cadherin (R&D Systems, AF648; 1:200); EEA1 (Santa Cruz, sc-137130; 1:100); EHD3 (LSBio, LS-C133741; 1:150); endomucin (Sigma-Aldrich, HPA005928; 1:100); eNOS (Abcam, ab76198; 1:200); ezrin (Cell Signaling Technology, 3145S; 1:300); FAM189A2 (Thermo Fisher Scientific, PA5-63414; 1:200); fibronectin (Abcam, ab2413; 1:200); FKBP51 (R&D Systems, AF4094-SP; 1:50); FXYD4 (Thermo Fisher Scientific, PA5-63570; 1:200); GFAP (Thermo Fisher Scientific, 14-9892-82; 1:200); glucocorticoid receptor (Cell Signaling Technology, 3660; 1:2,000); glutathione peroxidase 1 (R&D Systems, AF3798; 1:100); glutathione peroxidase 3 (R&D Systems, AF4199; 1:50); glycophorin A (R&D Systems, MAB1228-SP; 1:500); GRP78 (Proteintech, 11587-1-AP; 1:200); HB-EGF (R&D Systems, AF-259; 1:100); histone H3 (Cell Signaling Technology, 4499; 1:400); HMOX1 (Thermo Fisher Scientific, MA1-112; 1:200); HSD11B2 (R&D Systems, MAB8630-SP; 1:100); KIM-1 (R&D Systems, AF1750; 1:200); IBA1 (Thermo Fisher Scientific, MA5-27726; 1:500); IDH1 R132H (Dianova, DIA-H09, 1:200); IL-1RA (Abcam, ab124962; 1:200; specificity issues were raised by the provider after our experiments were completed. We kept it in the panel as none of our findings were affected and we did not perform any biological inferences on the basis of this antibody); iNOS (Thermo Fisher Scientific, MA5-41652; 1:200); integrin-α1 (R&D Systems, AF5676; 1:300); integrin-α3 (Proteintech, 66070-1-Ig; 1:2,000); integrin-β1 (Abcam, ab179471; 1:800); Ki-67 (Agilent, M724029-2; 1:200); laminin (Abcam; ab11575, 1:200); LAMP1 (Cell Signaling Technology, 9091; 1:300); LC3B (Cell Signaling Technology, 3868; 1:300); LEL-DyLight 649 conjugate (Vector Laboratories, DL-1178; 1:300); LTL biotinylated (Vector Laboratories, B-1325-2; 1:500); MCT1 (Thermo Fisher Scientific, MA5-18288; 1:300); MerTK (R&D Systems, AF591; 1:200); MPO (R&D Systems, MAB3174; 1:200); nephrin (Progen, GP-N2; 1:150); neurofilament (Agilent, IR607; 1:200); NOX4 (R&D Systems, MAB8158;, 1:300); NQO1 (Proteintech, 67240-1-Ig; 1:2500); OLIG2 (Bio SB, BSB 2561; 1:200); p62 (Cell Signaling Technology, 39749; 1:400); PCK1 (Proteintech, 66862-1-Ig; 1:400); PCNA (Abcam, ab29; 1:2,000); PDGFRβ (Cell Signaling Technology, 3169; 1:100); PDI (Cell Signaling Technology, 45596S; 1:400); periostin (R&D Systems, AF3548; 1:150); phospho-AMPKα (Cell Signaling Technology, 2535; 1:200); pJUN (Abcam, ab32385; 1:200); phospho-ERK1/2 (Cell Signaling Technology, 4370; 1:250); phospho-ezrin–radixin–moesin (Cell Signaling Technology, 3726; 1:200); phospho-GSK3β (Cell Signaling Technology, 9323; 1:100); phospho-histone H3 (Cell Signaling Technology, 9701; 1:200); phospho-JAK2 (Thermo Fisher Scientific, MA5-42424; 1:100); phospho-ribosomal protein S6 (Cell Signaling Technology, 4858S; 1:300); phospho-SMAD2 (Thermo Fisher Scientific, 44-244G; 1:200); phospho-SMAD3 (Thermo Fisher Scientific, PA5-104940; 1:200); phospho-STAT1 (Cell Signaling Technology, 9167S; 1:400); phospho-STAT3 (Abcam, ab76315; 1:200); PITX2 (R&D Systems, AF7388; 1:100); podocin (Sigma-Aldrich, P0372; 1:3,000); proteasome 20S LMP7 (Abcam, ab3329; 1:400); RAB5A (Cell Signaling Technology, 46449; 1:300); RAB7 (Abcam, ab137029; 1:200); RAP1GAP (Abcam, ab244259; 1:300); RCAS1 (Cell Signaling Technology, 12290; 1:200); sclerostin (Thermo Fisher Scientific, PA5-37943; 1:100); SIRT1 (Cell Signaling Technology, 8469; 1:200); SLC12A3 (Thermo Fisher Scientific, MA5-41643; 1:200); SOD1 (Proteintech, 67480-1-Ig; 1:400); SOD2 (Thermo Fisher Scientific, PA5-30604; 1:300); SRB1 (Abcam, ab217318; 1:300); STAT2 (R&D Systems, MAB16661; 1:200); survivin (Cell Signaling Technology, 2808; 1:300); talin 1 (Abcam, ab71333; 1:200); TRPC6 (Abcam, ab233413; 1:200); ubiquityl-histone H2B (Cell Signaling Technology, 5546T; 1:200); uromodulin (R&D Systems, AF5144; 1:300); villin 1 (Abcam, ab52102; 1:200); vimentin (Progen, GP53; 1:200); von Willebrand factor (Agilent, A008229-2; 1:200); WT1 (Agilent, IS05530-2; 1:200); and ZO-1 (Thermo Fisher Scientific, 61-7300; 1:250).

For mouse samples

ACE-2 (R&D Systems, AF933; 1:200); AIF (Cell Signaling Technology, 5318; 1:200); AKAP12 (Proteintech, 25199-1-AP; 1:600); ANXA3 (Sigma-Aldrich, HPA013398; 1:200); αSMA-FITC conjugate (Abcam, F3777; 1:800); aquaporin 2 (Alomone Labs, AQP-002; 1:400); calreticulin (Abcam, ab92516; 1:300); caspase 1 p20 (Thermo Fisher Scientific, PA5-99390; 1:200); CD3 (Abcam, ab1108; 1:200); CD4 (Abcam, ab183685; 1:200); CD41 (Thermo Fisher Scientific, PA5-79526; 1:500); CD42b (Abcam, ab227669; 1:100); CD44-Alexa Fluor 647 conjugate (BioLegend, 103018; 1:200); CD45 (Cell Signaling Technology, 70257; 1:200); FOS (Abcam, ab190289; 1:600); collagen I (Southern Biotech, 1310-01; 1:200); collagen IV (Abcam, ab6586; 1:200); cytochrome c (Abcam, ab110325; 1:200); DACH1 (Sigma-Aldrich, HPA012672; 1:200); E-cadherin (R&D Systems, AF648; 1:200); endomucin (Santa Cruz, sc-65495; 1:200); fibronectin (Abcam, ab2413; 1:200); histone H3 (Cell Signaling Technology, 4499; 1:400); IBA1 (Thermo Fisher Scientific, MA5-27726; 1:500); IL-1RA (Abcam, ab124962; 1:200; specificity issues were raised by the provider after our experiments were completed. We kept it in the panel as none of our findings were affected and we did not perform any biological inferences on the basis of this antibody); Ki-67 (Abcam, ab15580; 1:200); lamin B1 (Santa Cruz, sc-374015; 1:200); laminin (Abcam, ab11575; 1:200); LTL biotinylated (Vector Laboratories, B-1325-2; 1:500); nephrin (Progen, GP-N2; 1:150); PCNA (Abcam, ab29; 1:2,000); PDI (Cell Signaling Technology, 45596S; 1:400); phospho-ezrin–radixin–moesin (Cell Signaling Technology, 3726; 1:200); podocin (Sigma-Aldrich, P0372; 1:3,000); podoplanin (R&D Systems, AF3244-SP; 1:200); synaptopodin (Synaptic Systems, 163 004; 1:200); tyrosine hydroxylase (Cell Signaling Technology, 45648; 1:200); ubiquityl-histone H2B (Cell Signaling Technology; 5546T; 1:200); β-actin (Sigma-Aldrich; A5441, 1:1500); vimentin (Progen, GP53; 1:200); and von Willebrand factor (Agilent, A008229-2; 1:200).

Secondary antibodies and biotin-binding proteins

Secondary antibodies were diluted in a ratio ranging from 1:200 to 1:300. The following antibodies were used: goat anti-guinea pig IgG Alexa Fluor 488 (Thermo Fisher Scientific, A-11073); goat anti-guinea pig IgG Alexa Fluor 555 (Thermo Fisher Scientific, A-21435); donkey anti-mouse IgG Alexa Fluor 488 (Thermo Fisher Scientific, A-21202); donkey anti-mouse IgG Alexa Fluor 555 (Thermo Fisher Scientific, A-31570); donkey anti-mouse IgG Alexa Fluor 647 (Thermo Fisher Scientific, A-31571); donkey anti-rabbit IgG Alexa Fluor 488 (Thermo Fisher Scientific, A-21206); donkey anti-rabbit IgG Alexa Fluor 555 (Thermo Fisher Scientific; A-31572); donkey anti-rabbit IgG Alexa Fluor 647 (Thermo Fisher Scientific, A-31573); donkey anti-goat IgG Alexa Fluor 488 (Thermo Fisher Scientific, A-11055); donkey anti-goat IgG Alexa Fluor 555 (Thermo Fisher Scientific, A-21432); donkey anti-rat IgG Alexa Fluor 488 (Thermo Fisher Scientific, A-21208); donkey anti-rat IgG Alexa Fluor 555 (Thermo Fisher Scientific, A78945); donkey anti-sheep IgG Alexa Fluor 488 (Thermo Fisher Scientific, A-11015); donkey anti-sheep IgG Alexa Fluor 555 (Thermo Fisher Scientific, A-21436); streptavidin Alexa Fluor 488 (Thermo Fisher Scientific, S11223); and streptavidin Alexa Fluor 555 (Thermo Fisher Scientific, S21381).

Immunofluorescence in rat and human specimens

FFPE tissues were cut at a thickness of 2–3 μm, carefully affixed onto Superfrost Plus adhesion slides (Epredia, J1800AMNZ) and dried overnight at 37 °C. Subsequently, samples underwent sequential treatment involving triple immersion in xylene (10 min each) followed by an ethanol series (5 min each) consisting of three rounds of 100% ethanol, two rounds of 70% ethanol, one round of 50% ethanol and finally triple immersion in ddH₂O (5 min each). The immunostaining procedure mirrored the one used for PathoPlex samples but substituted 5% BSA with SuperBlock blocking buffer (Thermo Fisher Scientific, 37535) during the blocking step. Finally, after immunostaining, samples were mounted using ProLong Gold (Thermo Fisher Scientific, P36930).

Microscopy systems

For Fig. 2, images were acquired using a LSM 800 confocal microscope plus AiryScan (Zeiss, ZEN2.6) with the optimized ×63 objective (NA: 1.4). For Fig. 3, a Thunder Imager Live Cell and 3D assay (Leica Microsystems) fitted with a ×40 (NA: 1.10) or ×63 (NA: 1.40) objective was used to acquire images, which were processed using a computational clearing algorithm (Leica Microsystems)⁶⁰. The positional data of the imaged region for each sample were stored in Leica Application Suite X software (v.3.7.6, Leica Microsystems), which ensured consistent capture of the identical location for each cycle. For Fig. 4, a CellDiscoverer 7 with LSM 900 (Zeiss, ZEN 3.5 System) and AiryScan Multiplex fitted with a ×50 (NA: 1.20) objective and ×0.5 zoom was used to acquire images. Supplementary Table 8 summarizes the approximate microscopy times per experiment, considering image acquisition as the most important contributing factor. However, there are additional practical contributors, including chamber repositioning, movement delay between the ROI and data storage, which should be accounted for during implementation.

3D printing

Tinkercad (Autodesk; https://www.tinkercad.com) was used to create designs for the 3D-printed parts. The design of the headpiece was adjusted on the basis of a previously proposed design⁶¹. The BLTouch Cover Size Fixed was designed by louise_tguk on Thingiverse (https://www.thingiverse.com/thing:5013058). The chamber frame, the table for the chamber frame, the corner frame, the stage, the solution container, stands 1 and 2 for the solution container, the discard container, the base plate, the headpiece, the alignment guide, the BLTouch cover size fixed and the BLTouch cover box were printed using PLA filament 1.75 mm (Flashforge). The inner frame was printed using NinjaFlex TPU filament 1.75 mm (NinjaTek). The dewaxing container, the dewaxing container holder, the dewaxing carrier and the carrier handle were printed using PolyLite PETG filament 1.75 mm (Polymaker). An Ender-5 Plus printer (Creality) was used. The following settings were implemented in Ultimaker Cura (v.4.13.1; Ultimaker): nozzle size, 0.40 mm; layer height, 0.20 mm; wall thickness, 2.0 mm (PLA and PETG for containers), 1.2 mm (PLA and PETG for others) and 0.80 mm (TPU); top and bottom thickness, 2.0 mm (PLA and PETG for containers), 0.8 mm (PLA and PETG for others and TPU); nozzle temperature, 190 °C (PLA) and 235 °C (PETG and TPU); bed temperature, 69 °C (PLA), 75 °C (PETG) and 50 °C (TPU); fan speed, 100%; print speed, 60 mm s^–1 (PLA), 40 mm s^–1 (PETG) and 20 mm s^–1 (TPU); first-layer print speed, 20 mm s^–1 (PLA), 15 mm s^–1 (PETG) and 10 mm s^–1 (TPU); infill, 20% and zigzag; build-plate adhesion, brim. Masking tape was used to create an adhesive surface on the bed.

3D printer-based liquid-handling system

To prepare for the use of the liquid-handling system, several preliminary steps were required. These encompassed manual deparaffinization, antigen retrieval and mounting the Nexterion glass with sample sections onto the chamber frame. The deparaffinization process described above required the use of a dewaxing container, a container holder, a dewaxing carrier and carrier handle printed with PETG. After completion of dewaxing, the Nexterion glass with sections underwent antigen retrieval and washing procedures as outlined above. After washing, any excess wash buffer present at the edges of the glass was carefully removed. The Nexterion glass was then inserted into the bottom of the chamber frame, and its edges were securely sealed using silicone. It was crucial to allow the silicone to dry for a minimum of 15 min. To prevent the samples from drying out during this process, regular application of wash buffer to the samples was necessary while ensuring that the silicone did not become excessively wet. Once the silicone was completely dry, the inner frame was positioned in the frame and samples were covered with wash buffer. The following process used a liquid-handling system based on a 3D printer (Ender-5 Plus). For light shielding, the 3D printer was covered by a 3D Printer enclosure (Creality). The window on the front of the enclosure needed to be covered with an opaque material to shield the inside from light. The BLTouch built into the Ender-5 Plus was also partially shielded by attaching the BLTouch Cover Size Fixed. Our liquid-handling system was based on three different g-codes: (1) ‘BSA to Elu.gcode’, which automated the process from blocking with BSA during the initial cycle to elution and the pouring of imaging buffer; (2) ‘1st Ab to Img.gcode’, which automated the process of washing the imaging buffer, incubating the primary and secondary antibodies and pouring the imaging buffer; and (3) ‘Elu to Img.gcode’, which automated the process of washing the imaging buffer, elution and pouring the imaging buffer. The settings of the solution containers corresponding to each g-code are shown in Supplementary Table 9. The solution container stand consists of numbered sections ranging from 1 to 12, which are designated for container installation. Each solution-filled container was placed in the section with the corresponding number on the stand. The corner frames were inserted into the holes at the four corners of the table. Solution container stands 1 and 2 with solution containers, the discard container, the stage, the table and the chamber frame were placed on the base plate. Specifically, solution container stand 1 needed to be positioned at the front side of the base plate. To prepare for the operation of the 3D printer, the print bed was removed and autolevelling was disabled. Once each g-code was initiated, after calibrating the home position, the printer head moved slightly backward and the printer stage was lowered. The printer then paused for 60 s before resuming operation. During this pause, the headpiece and BLTouch cover box were attached to the printer head and BLTouch, respectively, and the base plate, complete with all the necessary components and solutions, was placed on the printer stage using the alignment guide. Once installation was complete, the enclosure was completely closed. The washing, staining and elution processes were then automatically performed by pushing down on the solution containers and table with a rod in the headpiece. The BSA to Elu.gcode, 1st Ab to Img.gcode and Elu to Img.gcode programs were completed in approximately 2 h 15 min, 3 h 10 min and 1 h 15 min, respectively. The dimensions of the chamber frame match those of ready-made plates used for imaging cultured cells. For an example of how this solution works, see Supplementary Video 1.

PEC cell line for migration assays

Primary PECs were thawed and cultured at 5% CO₂ and 37 °C in endothelial cell basal medium (ECBM; PromoCell, C-22210) and 20% FBS (Thermo Fisher Scientific; 10500064) until 70%–80% confluence. The maintenance culture was passaged three times a week by gentle trypsinization using trypsin-EDTA 0.05% (Thermo Fisher Scientific, 25300054).

Migration assays were performed using Culture-Insert 2 wells in μ-Dish 35 mm (Ibidi, 81176). Each well was seeded with 30,000 PECs in 100 μl ECBM without supplements and with 1% FBS and incubated overnight. The insert was then removed, which created a gap of 500 μm between cells. PECs were stimulated with either PDGF-BB (Thermo Fisher Scientific, 315-18-50UG, per manufacturer’s recommendations of 2.0 ng ml^–1) or with PDGF-BB and tanzisertib HCl (CC-930; Selleck Chemicals, S8490) or PDGF-BB and vehicle, in this case DMSO (Merck; D2650). All combinations were diluted using ECBM without supplements and with 1% FBS. Images were taken every 10 min for 23 h using the Personal Automated Lab Assistant (Leica Microsystems). Areas of migration were measured using Fiji.

Scratch assays were performed using glass-bottom FluoroDishes seeded with 50,000 cells in 2 ml ECBM with 1% FBS at 5% CO₂ and 37 °C. After 48 h, cells reached 90% confluence and were ready for the experiment, in which a sterile plastic 1,000 μl micropipette tip was used to scratch the cell monolayer to create a wound of around 1,000 μm. Next, the cell monolayer was gently washed with ECBM with 1% FBS to remove dead cell debris. To use the nucleus for tracking, PECs were stained with 80 nM Hoechst 33342 (Thermo Fisher Scientific) for 20 min at 5% CO₂ and 37 °C and washed once with Dulbecco’s PBS (Thermo Fisher Scientific, 14190094). Afterwards, 2 ml fresh ECBM with 1% FBS was added for image acquisition. Time-lapse imaging was performed using a Leica DMi8 M/C/A inverted microscope equipped with ×10 Plan Apo objective (Leica Microsystems). Images at both sides of the wound were acquired every 5 min with an ORCA-Flash4.0 digital camera (Hamamatsu Photonics) using MetaMorph (v.7.10.3.279) software (Molecular Devices). To visualize the wound, adjacent positions were stitched using the Stitching plugin from Fiji ImageJ. Tracking of the first 8 h of migration was performed with the TrackMate plugin from Fiji (v7.10.2) and custom-made scripts⁶². Mean square displacement was calculated using the CelltrackR package⁶³.

Transmission electron microscopy

For electron microscopy with immunogold labelling, kidneys were removed, cut into 2-mm-thick razor blade sections and immersion-fixed in freshly prepared 4% paraformaldehyde for 24 h at 4 °C. The samples were then resliced into 50-µm-thick sections using a vibratome. Vibratome sections were incubated with the primary antibody against TRPC6 (rabbit, final concentration 1:200). After washing and overnight incubation at 4 °C with the secondary antibody, goat anti-rabbit 1:100 (Nanoprobes), sections were silver enhanced with HQ silver (Nanoprobes) for 8 min in the dark at 4 °C, washed in 0.1 M phosphate buffer, treated with OsO₄ (0.5% for 45 min at room temperature) and stained with uranyl acetate (1% w/v in 70% v/v ethanol, 30 min at room temperature). After dehydration, sections were embedded in epoxy resin, Durcupan ACM (Sigma-Aldrich). Next, 50-nm ultrathin sections were cut using an UC6 ultramicrotome (Leica Microsystems) and analysed using an 80 kV Zeiss Leo 910 transmission electron microscope.

Imaging mass cytometry

Tissue sections were dewaxed in xylene and rehydrated, followed by staining using a standard protocol for immunohistochemistry according to the protocol by Fluidigm. Nuclei were labelled with iridium, and TRPC6 antibody (Abcam, ab105845) was coupled to 174Yb heavy metal. Data acquisition was performed on a Helios time-of-flight mass cytometer (CyTOF) coupled to a hyperion imaging system (Fluidigm). Areas for ablation were selected on the basis of haematoxylin and eosin staining performed on an adjacent slide. All data were collected using the commercial Fluidigm CyTOF software (v.01).

Pathological examination of crescentic nephritis

Images were evaluated by two expert pathologists in a blinded fashion to define disease states as either control, acute or crescentic phase. Next, multiple metrics were also calculated in a blinded fashion, including tubular injury score (0–3+), cell numbers per cross section, percentage of cellular crescents and percentage of tubular injury.

Bulk RNA-sequencing sample preparation and analysis

Glomeruli were isolated at day 4 after NTS treatment and in control groups after kidney perfusion with Dynabeads (Invitrogen), preserved in RNAlater and stored at −80 °C until processing. For preparation of nuclei, nuclei were extracted from the isolated glomeruli according to a modified protocol⁶⁴. The nucleus suspension was incubated on the magnet to remove magnetic beads used for the isolation of the glomeruli. Nuclei were mixed with RLT buffer (Qiagen) and frozen at −80 °C. Total nuclei RNA was extracted using RNeasy Micro kits (Qiagen) according to the manufacturer’s recommendations.

Bulk RNA-sequencing data were processed using our previously published open-source Snakemake⁶⁵ workflow for RNA-sequencing analysis with pytximport⁶⁶. In brief, raw FASTQ files provided by the sequencing facility were assessed for quality with FastQC⁶⁷, followed by trimming of adapter sequences and removal of low-quality reads with fastp⁶⁸. Next, processed sequences were selectively aligned to a reference gentrome based on Ensembl GRCh38 (release 112)⁶⁹ and transcript counts were quantified with Salmon⁷⁰. We used pytximport to estimate gene counts from transcript abundances with counts_from_abundance set to length_scaled_tpm. Differentially expressed genes were identified using PyDESeq2 (ref. ⁷¹) with log₂-transformed fold shrinkage. Genes were considered differentially expressed if their log₂-transformed fold change value was greater than 0.5 or lower than −0.5 and their false-discovery-rate-adjusted P value was less than 0.01. The results from PyDESeq2 were used to infer transcription factor activity based on the CollecTRI⁷² gene regulatory network reference with decoupler-py⁷³ univariate linear modelling.

Python library for spatial proteomics

Commonly used functions for the analysis of PathoPlex imaging data were combined into the spatiomic Python package. spatiomic comprises different submodules that facilitate data loading, image registration, image preprocessing, dimensionality reduction, spatial analyses, neighbourhood graph construction and clustering. In this section, we describe the architecture of this software package and the general functionality it includes. Parameter choices and detailed workflows are described in subsequent sections. Both the computational analysis library and the code for all analyses are available online through GitHub and Zenodo, as detailed in the Code and Data availability sections.

spatiomic comes with support for multiple common microscopy imaging formats and flexibly supports AnnData⁷⁴ objects, NumPy⁷⁵ arrays and pandas DataFrames⁷⁶. It uses cuml, cucim and cugraph from the RAPIDS ecosystem²¹ as well as cupy⁷⁷ for GPU-accelerated computations, which enables analyses to scale to billions of data points with affordable hardware with a time requirement of just minutes to hours depending on dataset size. The library was extensively unit-tested, typed and documented. Documentation, including a full example notebook detailing how to apply spatiomic analyses to PathoPlex data, is available at: https://spatiomic.org/.

To enable modular composition of analyses, spatiomic encompasses several submodules as described below.

Data submodule

The data.read class offers a method for parsing common microscopy imaging formats such as .tiff, .lif and .czi files through readlif, tifffile and aicspylibczi bindings. A random subsample of imaging data can be obtained through the data.subsample class. The submodule further contains functionality to subset multichannel images to a specified list of channels and to export data to AnnData⁷⁴ objects, which enables interoperability with the scverse²³.

Process submodule

The process submodule offers common functions for preprocessing imaging data. This includes the clip class for channel-wise histogram clipping to either absolute values or percentiles, the zscore class for channel-wise z scoring and the normalize class for the channel-wise scaling of intensity values to a given range. It also includes the register class, which exposes different methods for image registration and registration evaluation.

Dimension submodule

Dimensionality reduction is an important step in many single-cell and spatial omics analyses. spatiomic provides classes that help reduce the dimensionality of both the channel dimension and the data point dimension. The former is possible through the integration of the dimension.pca, dimension.tsne and dimension.umap classes, which internally rely on the Python packages scikit-learn⁷⁸, umap-learn⁷⁹ and cuml²¹. The latter is achieved by incorporating the dimension.som class, which enables GPU-accelerated training of SOMs thanks to an XPySOM²² integration. To train SOMs with the Pearson correlation as distance metric, we extended XPySOM with a CuPy-based function, which is available from GitHub (https://github.com/complextissue/xpysom).

Neighbour submodule

The neighbour submodule exposes classes that enable the creation of k-nearest and shared nearest neighbour graphs, which facilitates the construction of similarity-based graphs for graph clustering and distance-based neighbourhood graphs for spatial analysis.

Cluster submodule

Clustering algorithms enable the unsupervised identification of similar (protein co-expression) patterns, which facilitates automatic partitioning of complex signals into biologically meaningful groups. spatiomic.cluster includes classes for GPU-accelerated clustering with the Leiden graph clustering algorithm²⁴, k-means and hierarchical agglomerative clustering.

Spatial submodule

The spatial submodule incorporates functions for the explorative analysis of spatial distribution patterns in both immunofluorescence and clustered images, including global and local univariate and bivariate measures of spatial distribution based on PySAL⁸⁰. It further includes code for efficient join count quantification and spatial vicinity graph construction and interoperability tools for use together with PySAL⁸⁰, thereby ensuring compatibility with a wide range of spatial statistics applications⁸¹.

Tool submodule

The tool submodule contains utility functions for additional evaluation or analysis. It enables quantification of cluster abundance and identification of significantly differentially expressed clusters, calculation of mean immunofluorescence marker intensities per cluster and identification of cluster-defining immunofluorescence markers.

Plot submodule

spatiomic includes plotting functions based on matplotlib⁸² and seaborn⁸³ that facilitate visualizing common plots, for example, image registration metrics, SOM training-quality metrics, cluster projections, spatial adjacency graphs as well as cluster contributor histograms and volcano plots.

System and time requirements

spatiomic is designed to be flexible and adaptable to the scale of multiplexed imaging data and available computing resources. The only mandatory system requirement is Python (v.3.10) or higher. Although many individual functions run in seconds, and a complete exemplary workflow (excluding data download) can be completed in less than 3 min on a standard personal laptop, spatiomic substantially benefits from CUDA-enabled GPUs compatible with the RAPIDS ecosystem for larger datasets. At the scale of analysis presented in this paper, three specific steps took several hours to complete on a single GPU. First, image registration time scales linearly with the number of ROIs acquired and imaging cycles. Second, the time for SOM training scales linearly with the training sample size, the number of training iterations and the number of SOM nodes. Finally, the transfer of cluster labels from the trained SOM to all images in the dataset scales linearly with the number of acquired ROIs and the number of SOM nodes. Therefore, users are encouraged to evaluate and test parameter choices to optimize performance for their specific experimental design. For example, smaller, more homogeneous datasets may benefit from smaller training subsamples and SOM sizes, which will result in faster analysis. Moreover, users should select hardware that is appropriate for their analytical needs and desired turnaround time.

Image registration and autofluorescence correction

Image registration is the first step of every computational analysis of PathoPlex images. Given the cyclical nature of image acquisition and imperfect repositioning of standard microscopy systems, aligning signals from all imaging cycles to a common reference is a prerequisite for joint analyses.

We propose two different ways of aligning iteratively acquired immunofluorescence images: on the basis of either a nuclear or a pan-protein stain as a registration reference. Although nuclear markers are widely used for registration, a pan-protein marker enables the alignment of images and facilitates the delineation of the portion of the image covered by tissue, thereby lowering the computational burden for downstream processing and reducing the cluster annotation time by masking the empty background. To combat autofluorescence due to red blood cells (RBCs) in samples with abundant RBCs, we further used a RBC marker to specifically remove this source of noise. Detailed parameter choices and method descriptions are provided for the respective experimental datasets in the subsequent subsections. Differences between datasets, such as the use of NHS-E for foreground segmentation and the use of glycophorin A for erythrocyte segmentation in selected datasets, represent the iterative improvement of PathoPlex, with earlier experiments informing strategies for further scaling in subsequent experiments. In particular, the preprocessing steps described in this section refer only to the transformation of nuclear staining or NHS-E for registration purposes and do not indicate the transformation of marker intensities as used for further analysis, unless indicated otherwise.

Mouse CGN experiment

To align cycles from the CGN mouse dataset (Fig. 2), nuclear reference images (DAPI in cycles 1–17, DRAQ5 in cycles 18–42) were clipped image-wise to the range between the 1st and 99.9th intensity percentiles. Intensity values were z scored and normalized to the range 0–1. Gaussian blur with a sigma of 1 was then applied to the nuclear reference images to reduce noise and to improve comparability across cycles, followed by histogram matching of all subsequent images to the histogram of the reference image from the first cycle. These processed nuclear images were then aligned to the processed refence nuclear image obtained in the first imaging cycle. To this end, all intensity values below the 70th intensity percentile were set to 0 in both the reference nuclear staining and the staining to be aligned to increase the contrast between the nuclei (constituting less than 30% of the imaged pixels in all images) and the unstained background, followed by x and y offset detection between images with the phase_cross_correlation function of the Python package scikit-image⁸⁴, as included in spatiomic.process.register. All channels were projected onto the first cycle on the basis of the detected offset. Finally, the maximum offset detected between any two cycles was subtracted from all sides of all registered images, which resulted in square-sized and equally sized images.

Human advanced DKD experiment

For the study of specimens from individuals with advanced DKD and from individuals without diabetes (Fig. 3), all images were first corrected for camera distortion using the remap function of OpenCV⁸⁵ by linear interpolation based on a reference image of a micrometre-scale microgrid captured with the respective microscope objective. Histogram matching of the nuclear channel for each cycle to the reference nuclear channel of the first cycle was used to increase the similarity between nuclear images. Following that, the phase cross-correlation offset detection, as implemented in spatiomic.process.register.get_shift, was used to align all images to the first imaging cycle on the basis of the x and y offsets between the nuclear (DRAQ5) channels. For each field of view (FOV), the overlap across all cycles was calculated and the images were cropped to the respective overlap area.

Human early diabetes treatment experiment

This experiment (Fig. 4) leveraged NHS-E, a pan-protein stain, as a registration reference. To facilitate registration and further image processing, we applied local mean downscaling with a 2 × 2 pixel kernel to all imaging data. We then clipped the NHS-E histogram to the range between the 1st and 95th percentiles and normalized the NHS-E intensity to the range 0–1. After NHS-E preprocessing, Gaussian blurring with a sigma of 1, thresholding at the 70th percentile and phase correlation-based offset detection with an upsampling factor of 5 were applied to identify the x and y offset between NHS-E images from all cycles, using the first imaging cycle as the reference. After correcting for this offset, SIFT-based homography detection⁸⁶ was performed and images from subsequent cycles were warped to match the reference from the first cycle. Pre-registration and post-registration structural similarity index metrics (SSIMs) were evaluated for all registration pairs. Manual registration checks were performed for images with the following criteria: the SSIM decreased following registration; the SSIM after registration was lower than 0.3; any offset was greater than 300 pixels; or the SSIM after registration was lower than 0.8 and the SSIM improvement lower than 0.05. To account for RBC-derived autofluorescence, glycophorin A intensities across all images were clipped to the range between the 90th and 99.5th percentiles, and Otsu thresholding was used to binarize signals and to create RBC masks. Masks were smoothed through the removal of small objects <72 pixels, two iterations of binary dilation and hole filling. To remove signals from empty background without tissue, processed NHS-E signals were combined across all cycles by minimum projection (thus excluding areas with lifting at any point) for each image and Otsu thresholding was used to binarize signals. All images were then restricted to the area covered by the NHS-E mask but not included in the RBC mask. Autofluorescence correction was performed using secondary-antibody cycles acquired repeatedly throughout the experiment. For each primary-antibody imaging cycle, tissue autofluorescence was estimated by interpolating the signals from the nearest preceding and subsequent secondary-only cycles at the same wavelength, thereby accounting for minor fluctuations in autofluorescence. This estimated autofluorescence was then subtracted from the corresponding channels with lower intensity clipping at zero.

Other immunofluorescence samples

To achieve reliable and accurate results, we used an iterative registration framework called Elastix⁸⁷ and a Python wrapper package called PyElastix⁸⁸. To improve contrast and to mitigate the effects of varying signal strength across cycles, we applied histogram equalization to the nuclear channels. We also reduced the computational load by rescaling the images by a factor of 0.25. The first cycle was established as a reference, and subsequent cycles were aligned to it using normalized correlation as the optimization metric. The registration process involved iterative steps at six different resolution levels, with 1,000 iterations per level. We used the rigid Euler transform to account for x and y offsets as well as rotation. Further preprocessing was performed for brain samples with reduced contrast by truncating pixel intensity values to high percentiles. This approach focused on the sparsely available nuclei. All registrations were manually verified, and if individual registrations were unsuccessful, minor adjustments were made to the downsampling factor, spatial sample number or iteration number until satisfactory registration was achieved to maximize image inclusion.

Alignment of PAS staining to immunofluorescence images

Owing to differences in size and lens characteristics between the immunofluorescence and PAS images, additional processing steps were required to align the two modalities. First, lens-specific calibration matrices were used to remove optical distortions as described above for the human advanced DKD experiment. Next, the PAS images were rescaled to match the physical pixel resolution of the immunofluorescence images. As only greyscale images were used for registration, the PAS images were converted from RGB to greyscale using the OpenCV⁸⁵ Python library. For the greyscale version of the immunofluorescence image stack, three structural marker channels (DRAQ5, LTL and collagen IV) were selected and combined into a single channel using weighted addition. To enhance visual similarity between the greyscale PAS and immunofluorescence variants, the PAS image was inverted, and histogram equalization was applied to both the PAS and immunofluorescence images. To establish an initial alignment for registration of the immunofluorescence onto the PAS image, the PAS images were centrally cropped to match the smaller size of the immunofluorescence images (3,000 × 4,000 versus 2,048 × 2,048 pixels, respectively). Once the x and y offsets in this subregion were determined by the registration algorithm, the original multichannel immunofluorescence images were transformed accordingly. The full-size RGB PAS images were also cropped to the size of the transformed immunofluorescence images to facilitate overlaying of the registration results.

Suitability of registration reference markers

In a subsequent investigation, our goal was to evaluate the suitability of DAPI, DRAQ5 and NHS-E as registration reference labels. To achieve this goal, we acquired images of all three marker channels for each imaging spot in each cycle. For each image requiring registration, transformations based on each of the three markers were independently computed. The corresponding DRAQ5 channel of each image was adjusted to account for the detected transformation, and the structural similarity index measure was computed by comparing it to the DRAQ5 channel of the reference cycle. This approach ensured that the different registration references could be compared using DRAQ5 as a quality metric, which consistently provided reliable and high-contrast images. All registrations were manually reviewed, and in instances where registration was not successful, minor adjustments were made to the downsampling factor, the spatial sample number or the iteration number until a satisfactory registration was achieved. This process ensured that no image had to be discarded and that high-quality alignments were obtained.

Quality control of automated elution

Complete elution of bound antibodies between cycles is key to the cyclic acquisition of images. To quantify elution efficiency, secondary-antibody-only staining cycles were acquired throughout the experiments to establish a baseline autofluorescence profile. For each primary-antibody staining, the 50th and 99.95th intensity percentiles were quantified and compared with the secondary-antibody-only cycles acquired at the same wavelength using secondary antibodies directed against the host species of the primary antibody. For instances when at least one value was equal to or lower than any intensity from the secondary-antibody-only cycles (as may be the case for non-abundant markers, for example, markers for rare immune cells or phosphorylation states), the images were manually re-evaluated to include specific staining.

Generation and interpretation of protein co-expression clusters

Protein co-expression clusters are the standard output of PathoPlex analyses with spatiomic, which captures specific co-expression patterns of different proteins at the pixel level. These patterns are jointly identified for all imaging data from each respective PathoPlex experiment, which results in a consistent clustering across all samples and facilitates comparisons of spatial expression or co-expression of immunofluorescence markers. Identification of these pixel-level clusters is a multistep process that consists of weighted random subsampling of the data to ensure equal representations of all desired variables, signal preprocessing, training of a SOM to identify representative co-expression signals and finally similarity graph-based clustering. Once clusters were identified, their constituting signals were compared to all other clusters and projected in space to infer the biological processes they represent. Last, cluster abundance was quantified and compared across conditions to delineate regulated signals. We applied this overall concept to each experimental dataset with parameters specified as described in the subsequent subsections.

Weighted random subsampling

To limit bias due to different data sizes between samples and to reduce the computational burden of pixel-based clustering, a weighted subsample of approximately 5 million (10 million for the early human diabetes dataset) random pixel positions per imaging plate were sampled. For each imaging plate (n = 1 for all experiments, except for the human advanced DKD experiment, for which n = 2 imaging plates were used), the number of subsampled pixel positions was equally distributed across all disease states, and for each disease state, all samples were equally weighted. Finally, each FOV of a given sample was given the same weight in the subsample. For the early human diabetes treatment dataset, the subsample also considered equal numbers of pixels from periglomerular, glomerular and tubular images.

Histogram clipping and normalization

Immunofluorescence markers differ in intensity range and abundance, which therefore requires preprocessing steps to improve comparability. On the basis of weighted random subsamples, histogram clipping and range normalization classes contained in the spatiomic.process submodule were fitted. After channel-wise fitting of the classes on the random weighted subsample, all channels of all images of each respective dataset were transformed according to the established clipping limits and normalization and scaling settings. For mouse samples and the advanced DKD experiment, histogram clipping was performed based on the 50th (lower) and 99.7th (upper) percentiles for each respective marker in the subsample. Owing to the extensive marker panel for the larger early human diabetes treatment experiment, coupled with the increased sensitivity of the confocal microscopy system used for this dataset, histogram clipping was performed based on absolute intensity thresholds established by human expert annotation in a condition-blind fashion, with each marker evaluated on random patches extracted from random images, followed by normalization of the clipped images.

SOM fitting

Pixel-based clustering was used to isolate groups of similar immunofluorescence marker signals (clusters), which formed the basis of all downstream analyses. The first step towards this clustering is the fitting of SOMs to the dataset to reduce data-point dimensionality and to ensure computational feasibility of graph clustering as well as to improve representation of signals that are relatively rare in the training data. For each imaging plate (n = 2 for the human advanced diabetes experiment, n = 1 for all other datasets), a SOM was trained on the corresponding weighted random subsample. The SOM was initialized with a grid size of 500 × 500 nodes (400 × 400 for the early diabetes experiment) and used the Euclidean distance metric for the mouse CGN dataset and the advanced human diabetes dataset with the cosine similarity metric used for the early human diabetes experiment. The training process used spatiomic.dimension.som. For human samples, a final learning rate of 10⁻⁴ and a final Gaussian neighbourhood sigma of 10⁻³ (3 × 10⁻³ for the early human diabetes dataset) were used. For the mouse samples, the default settings for the learning rate and neighbourhood size were used. The training process was repeated for 50 iterations for all datasets.

Graph clustering and batch integration

On the basis of the representation of the signals contained in each experiment as provided by the respective SOM nodes, we used Leiden graph clustering to identify clusters of protein co-expression patterns and applied the clusters to all pixels from all images for each dataset. First, the similarity-based neighbourhood graph of SOM nodes was built using spatiomic.neighbor.knn_graph using cosine (early diabetes treatment experiment) or Euclidean distance (all other experiments). When image acquisition was performed using multiple plates (Fig. 3), the neighbourhood graph construction step was modified to use an adaptation of the batch-balanced k-nearest neighbours⁸⁹ algorithm, implemented in spatiomic.neighbor.knn_graph. A neighbour count of k = 40 for each respective plate was used. When only one imaging plate was used (all other experiments), a neighbour count of k = 40 (k = 50 for the early human diabetes experiment) was used without any batch integration. After graph construction, we used the Leiden²⁴ graph clustering algorithm to identify clusters of similar protein co-expression patterns. A Leiden resolution of 2.5 and an iteration count of 10,000 were used for the advanced DKD experiment, a resolution of 2.0 and an iteration count of 1,000 for the early human diabetes treatment experiment and a Leiden resolution of 1.0 and an iteration count of 10,000 for the mouse CGN experiment samples.

Cluster identity

To assign biological identities to protein co-expression clusters derived from PathoPlex output, we used a two-step approach that combined statistical overrepresentation with spatial distribution analysis and signal interpretation by human experts. Owing to computational constraints, statistical analyses were performed on subsampled data (used with the early diabetes experiment, with 500,000 pixels randomly selected from the weighted random subsample) or the representative SOM nodes (all other datasets). For each cluster, we calculated the mean normalized intensity and the log₂-transformed fold change of marker intensity (relative to all pixels or nodes not assigned to that cluster), which reflected the signal strength and cluster specificity, respectively. Significance was assessed using a two-tailed t-test with Benjamini–Hochberg (early diabetes dataset) or Holm–Šidák (all other datasets) correction for multiple testing (implemented in spatiomic.tool.get_stats). Dominant markers in each cluster (termed high contributors) were identified by ranking markers based on the product of their mean intensity and log₂-transformed fold change, and retaining only those with mean normalized intensity ≥ 0.2, adjusted P < 0.05 and log₂-transformed fold change ≥ 1. Intensity histograms further visualized marker distribution patterns in each cluster for all markers and were also considered for cluster interpretation. Individual clusters were spatially projected onto the corresponding immunofluorescence and—where available—PAS-stained images. Human experts then visually validated and refined cluster identities by assessing their spatial distribution, considering morphological structures and correlations with immunofluorescence signals and integrating contextual biological knowledge.

Selection of specific foreground clusters

As multiplexed imaging data from Figs. 2 and 3 did not include a pan-protein marker in their panel, no foreground segmentation was performed before pixel-level clustering, which resulted in multiple clusters that corresponded to empty background areas. To account for this limitation (largely circumvented by NHS-E and a larger antibody panel in Fig. 4), we restricted extended analyses (but not differential cluster abundance analyses) to specific foreground clusters only or we treated all background clusters as a single cluster. For Fig. 4, we similarly assessed the specificity of all clusters and quantified their relative frequency, focusing visualizations and biological assessment to clusters deemed to represent specific signal (excluding minor imaging artefacts, autofluorescence signals and imperfect foreground segmentation) and accounting for >0.1% of foreground pixels. The specificity of clusters, defined as the extent to which clusters correspond to clear biological processes and/or structures, both at the spatial and at the molecular level, was assessed individually by a panel of three experts. In controversial cases, the experts presented their arguments until a unanimous decision for exclusion was reached. If arguments for exclusion were not unanimously accepted, the cluster was not excluded.

Cluster abundance and differential abundance analysis

Cluster abundances were analysed by quantifying the number of pixels assigned to each cluster for each field of view in a dataset using spatiomic.tool.count_clusters. The abundances were normalized to a range of 0–1 to facilitate comparison and interpretation across FOVs. For the mouse CGN and the advanced diabetes experiment, these normalized abundances were further aggregated to obtain the mean cluster abundance per field of view for each mouse or patient. Differential abundance analysis was performed by quantifying the log₂-transformed fold change in cluster abundance between different experimental and clinical conditions by performing a two-sided t-test. To account for the higher number of unique clusters and larger sample sizes, P values for the advanced DKD and the early diabetes treatment experiments were corrected for multiple testing with Benjamini–Hochberg false discovery rate adjustment. Statistical testing used the spatiomic.tool.get_stats function. Results were visualized with spatiomic.plot.volcano for all datasets. For visualization purposes, the volcano plots were restricted to feature clusters determined to represent specific foreground signal. Moreover, to assess the impact of batch integration, quality control was performed by comparing differences in log₂-transformed fold changes in inter-group cluster abundances between the imaging plates (Supplementary Fig. 12).

Extended analyses based on pixel-level clusters

Based on our foundational pixel-level protein co-expression clusters, we used multiple downstream applications to connect the output from PathoPlex with community resources and public knowledge bases, showcasing how pixel-level data can be aggregated at different levels to derive information across multiple biological scales.

Biclustering

UnPaSt⁴² is a biclustering method initially developed for unsupervised patient stratification based on omics data. UnPaSt identifies differentially expressed biclusters in a two-dimensional matrix with samples (for example, images or patients) in columns and features (for example, clusters) in rows. A differentially expressed bicluster is a submatrix consisting of samples and features such that these features are overexpressed or underexpressed in these samples compared with all other samples in the input data matrix. We applied UnPaSt with binarization P value threshold of P = 0.05, direction = ‘both’ (to identify biclusters consisting of both upregulated and downregulated clusters) and all other parameters set to default to image-level and patient-level cluster intensities. As UnPaSt is not deterministic, consensus biclusters were built on the basis of results of ten independent runs.

Druggability profiling

To evaluate the druggability of the molecular signature of advanced DKD and to identify potential therapeutics, we combined drug–protein interaction data from the CTD⁴⁴ with STRING⁹⁰ protein–protein interaction data. Initially, we manually identified AIFM1, TRPC6, CALR, HSPA5, ITGB1 and CTNNB1 as possible targets involved in the altered signalling cascades revealed by PathoPlex based on the significantly differentially expressed clusters and their molecular composition. Next, we used STRING data for Homo sapiens to link every potential target to a broader network of interacting proteins. We only considered direct interactions with a confidence score exceeding 0.75. In a subsequent step, CTD chemical–gene interaction data were used to extract possible therapeutics that affect proteins in our target networks. This search was limited to compound–protein interactions described in mice (Mus musculus), rats (Rattus norvegicus) or humans. Only interactions that did not involve co-treatment and did not affect the reaction of another externally administered compound were included. Our goal was to determine the impact of existing pharmacological treatments on the extended protein networks and to explore the potential for repurposing drugs authorized for other indications; therefore, we further filtered the results, preserving only entries for which the chemical name had a matching entry in the European Medicines Agency’s list of authorized agents (date of consultation: 23 August 2023). In a second step, to further assess the impact of current antidiabetic drugs, results from all protein networks were combined and filtered for compounds containing ‘glutid’, ‘gliflozin’, ‘glitazon’, ‘gliptin’, ‘metformin’, ‘pril’ or ‘sartan’.

Pixel cluster-assisted cell-level metaclustering

Although pixel clusters provide valuable subcellular and extracellular information, they can also be used to inform existing cell-level clustering workflows. To quantify the co-occurrence of pixel-level clusters within cell-level metaclusters, we first applied the Cellpose³⁹ segmentation model using the parameters model_type = “nuclei” and diameter = 30 to the pre-processed DRAQ5 channel of each image. Statistical testing of nucleus counts per image was performed using statannotations⁹¹ with a two-sided nonparametric Mann–Whitney U-test, comparing diabetes and control samples. Centroids of all identified nuclei were then expanded radially by the smaller of either 5 µm or 50% of the shortest Delaunay triangulation edge length to approximate cell areas. Within these estimated areas, the relative abundances of pixel-level clusters were calculated to produce feature vectors for each cell. To identify cell-level metaclusters, a k-nearest neighbour graph was constructed using spatiomic.neighbor.knn_graph with neighbor_count = 50 followed by Leiden²⁴ clustering with a resolution of 1.0. Condition-specificity of meta-clusters was confirmed by averaging their relative abundances across images at the patient level. Within each cell-level metacluster, the fraction of pixels corresponding to each pixel-level cluster was quantified, which provided a quantification of the cell-level co-occurrence patterns of pixel clusters. Differences in normalized metacluster abundance at the patient-level between conditions were visualized through a two-component principal component analysis using spatiomic.dimension.pca with default parameters and significance was established through a Mann–Whitney U-test with Bonferroni correction.

Pseudotime analysis

PILOT⁴¹ is a previously published, multiscale, unsupervised method that uses optimal transport to compute the distance between data points and infer a disease trajectory. First, we used the normalized abundance of the clusters as the cluster proportions for each FOV from all samples. Next, we computed the distances between the clusters as a cost matrix. In the subsequent step, we used the cluster proportions and the cost matrix to compute the Wasserstein distance⁹² (W1) between the data points. Then, we obtained the trajectory of the disease by applying the diffusion map⁹³ to the distance matrix of the samples. Finally, we used the assigned pseudotimes of the data points to reveal the changes in the proteins or clusters. In summary, PILOT uses stepwise nonlinear models to determine significantly changing proteins or clusters across the disease trajectory.

Cluster join counts analysis

The immediate spatial neighbourhood analysis was based on join counts between unique clusters that were used to create an adjacency graph using spatiomic.spatial.vicinity_composition. First, for each pixel position, the eight surrounding pixel positions (Chebyshev distance of 1 pixel; that is, first-order queen neighbourhood) were examined to count the instances of nearby clusters. This process was applied to all pixels in all images using a vectorized approach. The cluster counts were then aggregated at the dataset level for images of the same disease condition. Next, connections between identical clusters were discarded, and the remaining connections were normalized to a range of 0–1. Non-foreground clusters were discarded. The resulting adjacency matrix was used to construct a directed graph for each condition, with the graph representing the relationships between clusters. Although the entire adjacency matrix was quantified and evaluated, a neighbourhood cluster abundance of 7.5% was established as minimum value to be included in the graph plot for visualization purposes, focusing the visualization on the most common adjacencies. The graph layout was calculated using the software packages Graphviz⁹⁴ and NetworkX⁹⁵.

Condition-specific structural patterns

Condition-specific structural patterns were identified using MISTy⁴⁰ (mistyR v.1.6.1). To that end, image data (112 images from advanced DKD samples and 310 from control samples) were aggregated at two different resolutions by summing the cluster counts in bins with a side length of 10 μm (62 × 62 pixels). Clusters capturing empty background, unstained tissue parts and nonspecific signals were collapsed into a single background cluster. To account for truncated bins at the edges of slides, we transformed the counts into proportions. We used these cluster proportions per bin as an intrinsic representation of the structure in a bin (MISTy intraview). To capture the broader tissue structure, we constructed the paraview by summing up the cluster proportions of the 20 nearest neighbours using family = “constant”, l = 20. To construct the paraview for the high-resolution aggregation, we computed the weighted sum of the cluster proportions of the 80 nearest neighbours using a Gaussian kernel with a bandwidth of 2.5 μm (corresponding to 15 pixels) (family = “gaussian”, l = 2.5, nn = 80). With these view compositions per aggregation, a MISTy model was independently trained for each sample. The MISTy models identified significant structural patterns in the different spatial contexts by associating the proportion of pixels belonging to each cluster in each spatial context to the target proportions in the intraview. MISTy can learn both simple linear relationships (for example, cluster X has a higher proportion if cluster Y has a lower proportion) and complex nonlinear relationships. By combining the predictions from the intraview and paraview for each cluster, MISTy enabled us to disentangle whether the prediction for a given cluster improves, and to what extent, when taking different spatial contexts into account. To compare the MISTy importance scores between conditions, we first computed the mean results per sample due to differing numbers of imaged FOVs. We then aggregated the MISTy results per patient and finally per condition (advanced DKD and controls). For each level of aggregation, group and view, we generated a graph representing the inferred relationships between clusters. In each graph, the nodes represent the clusters and the edges between the clusters were weighted by the importance scores inferred by the MISTy model (thresholded to conserve only significant relationships with importance > 1.0). The graph layout was calculated using Graphviz⁹⁴ and NetworkX⁹⁵.

Single-cell and single-nucleus RNA sequencing

To contrast our findings at the protein-level with transcriptomic data and to establish a bridge to proposed treatments of DKD, we leveraged two public RNA-sequencing datasets with pharmacological intervention or different treatment data, covering rodent and human samples.

Processing single-cell RNA-sequencing data

Data were downloaded from the Gene Expression Omnibus (accession: GSE220939). Files from individual patients were converted to AnnData⁷⁴ objects, and information on the diabetes and SGLT2i treatment status of each patient was added. To remove ambient RNA contamination, Cellbender⁹⁶ (v.0.3) training was performed on each sample individually with a training fraction of 0.5 for 100 epochs. Observations with a Cellbender cell probability ≤ 0.5 were discarded. Next, all observations were combined into a single AnnData object, quality-control metrics were quantified and cell clustering was performed using the Python package scanpy⁹⁷. As part of this step, mitochondrial, ribosomal and haemoglobin genes were identified, and quality control metrics were calculated using the calculate_qc_metrics function. Barcodes containing ≥50% mitochondrial RNA, ≥20% rRNA or ≥5% haemoglobin genes were discarded. Barcodes with reads for less than 500 or more than 5,000 genes were discarded to correct for doublets, as were genes detected in at most 4 observations. The total number of reads per cell was normalized and the counts were log1p transformed. To reduce dimensionality, the first 50 principal components were computed using scanpy’s pca method with the arpack singular value decomposition solver. As inter-individual batch effects were present, Harmony⁹⁸ data integration was performed with standard parameters, using the harmonypy⁹⁹ implementation available through the external module of scanpy. On the basis of the adjusted principal components, the 50 nearest neighbours of each cell were identified using the neighbours function of scanpy⁹⁷, which was configured to use the Pearson correlation similarity metric. Graph clustering was performed using the Leiden algorithm⁷² with a resolution of 0.5. On the basis of cell-identity genes (adapted from a previous study⁴⁶ and PanglaoDb¹⁰⁰), clusters were manually annotated to represent different cell types. Cell types and transcripts related to the pathways altered at the protein level were chosen, and gene counts were visualized using the dotplot function of scanpy⁹⁷. Statistics per cell type of interest were calculated on the depth-normalized gene counts using the get_stats function of spatiomic, configured to revert log1p transformation. Nonparametric independent Wilcoxon rank-sum testing was performed using the ranksum function of scipy¹⁰¹ with Benjamini–Hochberg correction.

Processing single-nucleus RNA-sequencing data

Preprocessed data were downloaded from the Gene Expression Omnibus (accession: GSE209821). Gene counts were read into R using Seurat and converted to AnnData⁷⁴ with SeuratData and SeuratDisk. In the next step, the AnnData object was read using the scanpy⁹⁷ Python library, and provided metadata were combined with the gene counts. The provided counts had already been filtered for quality-control metrics and cluster labels from a previous study were provided⁴⁷, thus, only depth normalization to the median read depth was performed. Observations were further filtered to only retain proximal tubule cells from controls, untreated mice or mice treated with soluble guanylate cyclase activators. Differential expression statistics of genes corresponding to defining markers of clusters upregulated in the proximal tubular compartment in the advanced DKD PathoPlex experiment were calculated based on depth-normalized gene counts using the spatiomic.tool.get_stats function, configured to run nonparametric independent Wilcoxon rank-sum testing using the ranksum function of scipy¹⁰¹ with Holm–Šidák correction.

Signal intensity of TRPC6 and AIFM1 in Btbr
^ob/ob mice

To further evaluate the expression of TRPC6 and AIFM1 in the proximal tubule of kidneys with metabolic damage, we assessed their expression in Btbr^ob/ob mice and control mice at 12 and 24 weeks of age. Immunofluorescence was performed on 5 FOVs per sample, with n = 3 samples per group and time point, which produced a total of n = 12 samples and 60 FOVs. The proximal tubule area was segmented by thresholding the LTL channel to an intensity greater than 1,500, followed by filling holes smaller than 1,000 pixels and removing objects smaller than 20,000 pixels. After proximal tubule segmentation, a random subsample of 50,000 pixel positions was selected from the proximal tubule area of each FOV. These subsamples were combined for each unique condition and time point. Changes in TRPC6 and AIFM1 expression were visualized through histograms of marker intensities, grouped by time point and condition. Statistical analysis was performed using the statsannotations⁹¹ Python package, for which the subsampled intensities served as input, with the corresponding condition as the label for each intensity. An independent t-test based on the statannotations⁹¹ Python package was used to assess statistical differences in TRPC6 expression between conditions at each time point.

General statistical analysis

Statistical analyses, including cluster composition, differential cluster abundance and differential gene expression, were performed using spatiomic.tool.get_stats and statannotations. These tools internally rely on scipy¹⁰¹ and statsmodels¹⁰² and were used to perform either Benjamini-Hochberg or Holm–Šidák-corrected two-tailed t-tests (for cluster composition and differential abundance analysis), Mann–Whitney U-tests (for differences in nucleus counts and cell-level metacluster abundance) or nonparametric Wilcoxon rank-sum tests (for differential gene expression from single-cell and single-nucleus RNA-sequencing data). Bulk RNA-sequencing differential gene expression analysis was performed using PyDESeq2, combining single-factor analysis using Wald tests with log₂-transformed fold change shrinkage with approximate posterior estimation generalized linear models¹⁰³. Detailed descriptions of these tests, including information on input data and additional filters, are provided throughout the Methods. All other statistical analyses, including the quantification of changes in cell migration, clinical parameters (for example, proteinuria or eGFR), sample-specific variables (such as age), PEC activation, injury patterns in CGN and principal component analysis, were performed using GraphPad Prism (v.9). Violin plots report median and interquartile values. Significance was evaluated using the unpaired t-tests with Welch’s correction comparing two continuous variables, a paired t-test for before and after settings, and the Brown–Forsythe, Welch ANOVA and Dunnett’s tests when comparing three continuous variables. Correlation analyses were performed using Spearman rank coefficients. Principal components were selected on the basis of the percentage of total explained variance using normalized cluster abundances at the image-level or the patient-level as input. Statistical significance was defined as P < 0.05 for all analyses, with a threshold of P < 0.01 applied for specific cases as outlined throughout the methods.

Reporting summary

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