Human tissue samples
Written informed consent was obtained from adult patients undergoing elective or emergency IBD surgery in accordance with NHS Research Ethics Committee (REC) approvals (TIP-18/WM/0237; GI Biobank-16/YH/0247; IBD Biobank-09/H1204/30). Healthy control samples were sourced from patients undergoing elective colorectal cancer surgery (uninvolved tissue) or stoma reversal. Patient metadata are summarized in Supplementary Table 1.
Archived diagnostic formalin-fixed paraffin-embedded (FFPE) blocks for ST and validation were obtained from Oxford University Hospitals NHS Foundation Trust via the Oxford Radciffe Biobank (ORB; REC:19/SC/0173) and the Oxford Centre for Histopathology Research. Additional FFPE blocks were obtained from the Friedrich-Alexander University Erlangen-Nuremberg (ethics: 23–131 bp).
Sample collection, handling and processing
Full-thickness resection tissue samples, surplus to clinical needs, were obtained by the operating surgeon, who identified involved and uninvolved tissue macroscopically. Pathology was confirmed through clinical histopathology reports. Samples were stored in Dulbecco’s modified Eagle’s medium (DMEM, Gibco) supplemented with 10% fetal calf serum (FCS; Sigma-Aldrich), 100 U ml−1 penicillin, and 100 U ml−1 streptomycin (Sigma-Aldrich) and transported to the laboratory on ice.
Tissue samples were rinsed with phosphate-buffered saline (PBS, Gibco) and examined for anatomical landmarks. Blunt dissection was performed on ice to isolate full-thickness (FT), submucosa (SM) and mucosa (MUC) layers. Following dissection, samples were immediately preserved in CryoStor CS10 (Sigma-Aldrich) and stored in liquid nitrogen for subsequent digestion following the manufacturer’s protocol.
For FFPE sections, tissues were oriented in cassettes, fixed in 10% buffered formaldehyde (VWR International) for 48 h, and transferred to 70% ethanol before being processed through a graded alcohol series for paraffin embedding. FFPE blocks were stored at −20 °C for long-term use until required for histological staining, immunohistochemistry or ST experiments.
Cell isolation, staining and sorting
Full thickness tissue samples were dissociated into single-cell suspensions for scRNA-seq using optimized protocols for epithelial and lamina propria (LP) fractions.
Epithelial cell enrichment
Epithelial cell dissociation was performed by adapting previously described protocols12,53. In brief, cryopreserved FT or MUC tissue samples were thawed for a minute in a 37 °C water bath then washed in 15 ml of supplemented DMEM before centrifugation. Tissue samples were cut into small fragments followed by further centrifugation. Samples were incubated in pre-warmed chelation medium (HBSS (Lonza) supplemented with penicillin and streptomycin, HEPES, FCS, 5 mM EDTA (Invitrogen) and 2 mM dithiothreitol (Thermo Fisher)) for 20 min, vortexed periodically, and then centrifuged to isolate crypt-containing supernatants. The chelation protocol was completed a total of two times to increase yield of isolated crypts. Recovered crypts were digested with TrypLE Express (Gibco) and DNase (Sigma) for 30 min to achieve a single-cell suspension, filtered through 70-µm and 40-µm filters to remove debris, and counted for viability. Between 0.5 and 1 million highly viable cells per sample were taken forward for cell staining. Cells were incubated with TruStain FcX (BioLegend) for 10 min at 4 °C, to block non-specific binding. Cells were then stained with a cellular indexing of transcriptomes and epitopes by sequencing (CITE-seq) antibody mix (CD4, CD8, CD45RO, CD45RA, PD1, CD103, CD56 and CD3, each at 1:100) (BioLegend), to identify and exclude any contaminating immune subsets or intra-epithelial lymphocytes and also stained with hashtag-oligo (HTO) antibodies at 0.75 µl per 1 × 106 cells per 100 µl staining volume (BioLegend), for 30 min at 4 °C. After staining, cells were immediately spun down, and resuspended at the desired pooling concentrations prior to loading on the 10X Chromium platform.
Immune (CD45+) and stromal (CD45−EPCAM−) cell enrichment
LP fractions were isolated using an optimized Miltenyi Biotech LP Dissociation Kit protocol tailored for FT ileal tissue. Cryopreserved samples were rapidly thawed and dissected into small fragments and resuspended in buffer L (DMEM (high glucose), FCS and HEPES). Tissue fragments underwent digestion at 37 °C for 2 h, with periodic mechanical dissociation using a blunt needle and syringe. Following digestion, single-cell suspensions were filtered through a 100-µm filter, washed with DMEM followed by centrifugation and counted to ensure high viability prior to further processing.
For staining and sorting, single-cell suspensions were incubated with antibody panels targeting markers to identify immune (CD45+ (1:100)) and stromal (CD45− EPCAM− (1:50)) populations following TruStain FcX (BioLegend) blockage. TruStain FcX was added at 5 µl per 1 × 106 cells in 100 µl staining volume. Additional CITE-seq and hashtag antibodies were used, as above, to enable multiplexing for scRNA-seq. DAPI stain (1:1,000) just prior to FACS sorting for live-dead differentiation. Flow cytometry (BD LSR II) and fluorescence-activated cell sorting (BD FACSAria III and BD FACSAria Fusion (BD FACS Diva Software), BD Biosciences) were employed to enrich cell populations. The accuracy of sorting gates was confirmed using compensation controls generated using single colour controls and fluorescence minus one samples. Live cells sorted into stromal (CD45−EPCAM−) and immune (CD45+) subsets were collected in Eppendorfs containing FCS and immediately spun down, counted, resuspended and pooled prior to loading onto the 10X Chromium platform.
LeviCell-based stromal cell enrichment
The LeviCell platform (LevitasBio) offered an alternative method for stromal enrichment. In brief, FT tissue samples underwent both epithelial crypt chelation and LP digestion. LP-derived single-cell suspensions were then stained with respective hashtag-oligo antibodies, counted and pooled in desired one-to-one ratio. Levitation buffer and CD45 depletion beads were added to the stained LP-derived single-cell suspensions per the manufacturer’s instructions (1004001). The suspension was loaded into the LeviCell-1.0 system, which employs paramagnetic levitation to separate cells based on viability. Live stromal cells (crypt-depleted (CD45−)) were enriched from the top output well, and dead and CD45+ cells were removed via the bottom output well. The isolated stromal fraction was directly processed for scRNA-seq.
Droplet-based scRNA-seq
Droplet-based scRNA-seq was undertaken using the 10X Chromium Single Cell platform in accordance with the manufacturer’s instructions (10X Genomics, 5’ v1.1 chemistry, CG000208, Rev F).
Approximately 30,000–40,000 cells were loaded per pool. The workflow involved encapsulating individual cells within gel beads in emulsion (GEMs), barcoding, reverse transcription of complementary DNA (cDNA) within GEMs, cDNA clean-up, amplification, and library construction for gene expression (GEX) and antibody-derived tags (ADT) data.
To enhance sequencing quality and depth, the workflow was modified by incorporating Jumpcode CRISPRclean Single Cell RNA Boost Kit (KIT1018 v1.0, Jumpcode Genomics). Post-ligation clean-up of GEX libraries included CRISPR RNA depletion using the Cas9–single guide RNA ribonucleoprotein complex. This step targeted unaligned reads, highly expressed ribosomal and mitochondrial genes, and non-variable genes for depletion. Following this, library construction proceeded in line with the manufacturer’s protocol.
Library quality and integrity were assessed using Agilent Bioanalyzer TapeStation to ensure compliance with sequencing requirements. Final libraries were pooled at a 4 nM concentration and sequenced on an Illumina NextSeq 500/550 (High Output v2.5, 150 cycles) or outsourced to a NovaSeq PE150. Sequencing depth and run parameters followed the specifications provided by the manufacturer (CG000208, Rev F).
Spatial transcriptomics
Five to ten 5-µm scrolls were sectioned from each FFPE tissue block and stored at −80 °C for RNA quality assessment. Total RNA was extracted using the Qiagen RNeasy FFPE Kit (73504) following the manufacturer’s protocol. RNA integrity was assessed using the RNA Pico Assay (Agilent), which quantified DV200 values. Samples with DV200 values above 50% were considered suitable for ST.
All FFPE blocks were stained with H&E to evaluate tissue morphology and pathology representation. Only samples meeting both RNA quality and tissue morphology criteria were selected for spatial gene expression analysis.
Visium FFPE protocol
The sample preparation workflow for 10X Genomics Visium Spatial Gene Expression (CG000408 and CG000409 rev. D) was followed, with specific optimizations. FFPE blocks were stored at −20 °C, slide drying time was extended to 48 h at room temperature, and deparaffinization and alcoholic rehydration steps were modified to minimize tissue detachment.
For the ST experiment, 5-μm sections were placed onto Visium FFPE spatial GEX slides, which feature 4 capture areas (42.25 mm2 each) containing 5,000 gene expression spots per area. Tissue sections were deparaffinized, stained with H&E (imaged using a Zeiss AxioScan at 10× and 20× magnifications), and de-crosslinked. Tissue digestion, gene-specific probe hybridization and cDNA synthesis were performed directly on the slide, incorporating spatial barcoding, per the Visium protocol (CG000407, rev. D). cDNA quantification via qPCR (KAPA SYBR FAST) informed the number of DNA amplification cycles. Libraries were then constructed following the manufacturer’s instructions. Library quality was verified using the Bioanalyzer High Sensitivity DNA Kit (Agilent) and concentrations were measured using the Qubit dsDNA HS Assay Kit (ThermoFisher). Libraries were pooled at 4 nM, with sample proportions adjusted based on the ratio of tissue section area to capture area.
Sequencing was performed on an Illumina NextSeq 500 platform (High Output v2.5, 150 cycles) with the following run parameters: read 1 (28 cycles), i7 index (10 cycles), i5 index (10 cycles), and read 2S (50 cycles).
Xenium FFPE in situ protocol
Slides for Xenium in situ were prepared following the manufacturer’s instructions with modifications. In brief, FFPE tissue blocks were sectioned at 5-μm thickness and placed on Xenium slides with a capture area of 235 mm2 (CG000578, rev. C). After drying, sections were deparaffinized and de-crosslinked (CG000580, rev. C) to facilitate probe hybridization. Custom-designed Xenium probes were hybridized to the tissue overnight, followed by probe ligation and annealing of rolling circle amplification primers. Sequential washing steps were performed to minimize autofluorescence and nuclei were stained with DAPI for imaging (CG000582, rev. E). Imaging was conducted using the automated Xenium Analyzer (CG000584, rev. E) which identified transcript counts through unique optical signatures for each gene. Transcripts with quality scores (Q-score) exceeding 20 were retained for downstream analysis. Cell boundaries were initially defined using DAPI images (see ‘Computational analysis’ for cell segmentation details). Xenium Prime analysis was carried out as above, with the addition of cell boundary staining per the user guide (CG000749, rev. B).
After the Xenium run, slides underwent H&E staining (CG000613, rev. A) to enable direct alignment of gene expression data with histological information from the same tissue section. H&E slides were imaged at 10× and 20× magnifications using the Zeiss AxioScan Z1.
Histopathological characterization and validation
Histopathological analysis and annotation
For ST analysis, corresponding digitized H&E images were reviewed. Two gastrointestinal pathologists (M.V. and E.F.) independently examined the slides with limited prior knowledge of patient clinical status. Manual annotations captured cellular details, pathology and tissue architecture. Both pathologists reviewed the same slides, and minimal inter-observer variability was observed, ensuring robust assessments. A comprehensive review of the formal histopathology report, including macro- and microscopic details, supplemented the findings.
Haematoxylin and eosin staining
Paraffin-embedded tissue sections (5 µm) were placed on Superfrost slides (Avantor) and dried. Slides were deparaffinized and rehydrated through a graded ethanol series. Sections were stained with modified Harris haematoxylin (Sigma) and eosin Y solution (Sigma), with staining durations adjusted for optimal colour development. After air-drying, slides were mounted using Leica mounting medium and coverslips. Stained sections were digitized at 10× and 20× magnifications using the Zeiss AxioScan Z1.
Immunohistochemistry
Paraffin-embedded tissue sections (5 µm) were mounted on Superfrost slides and pre-heated at 60 °C for 1 h. Slides were deparaffinized in 100% Histoclear and rehydrated through graded ethanol. Heat-induced antigen retrieval was performed in pH 6 or pH 9 buffer (EnVision FLEX, Dako) using a Decloaking Chamber. Endogenous peroxidase activity was blocked with H2O2, and non-specific binding reduced using 2.5% horse serum.
Slides were incubated with primary antibodies (F3, pH 6 1:250; TWIST1, pH 6 1:800; CD45, pH 6 1:50) for 1 h at room temperature or overnight at 4 °C. After washing with Flex Buffer (Dako), slides were treated with horseradish peroxidase-conjugated secondary antibodies (Dako) and washed again. Peroxidase activity was visualized using DAB substrate (Dako), and slides were counterstained with haematoxylin for nuclear detection. Following washing, slides were dried and mounted with Leica Mounting Medium. High-resolution images (20× magnification) were captured using the Zeiss AxioScan Z1.
Collagen fibre staining and imaging
FFPE sections (4 µm) were deparaffinized through a standard xylene and ethanol series to water. They were incubated with Picrosirius Red (HB6179) for 60 min at room temperature. Once excess was blotted off, slides were washed twice in 200 ml 0.5% acetic acid solution for 1 min each, followed by 100% ethanol washes for 4 min each. After dipping in Histoclear, slides were covered with micromount medium and coverslip. Imaging was performed at 40× magnification on a Zeiss Axioscanner (brightfield imaging). Polarized light imaging was performed on an Evident VS200 slide scanner (polarized light illumination), using a 40× 0.95 NA air objective lens, with an iDS VS-264C colour camera. An automated stage was used to tile and stitch images of tissue sections using the VS200 ASW software.
Quantitative PCR validation
RNA from FFPE tissue blocks of fistulating CD patients and controls (CD, non-CD and healthy individuals) was reverse-transcribed into cDNA using the Reverse Transcription Master Mix (1006300, Standard Bio). TaqMan assays (Thermo Fisher) were pre-amplified, diluted, and used for qPCR.
The Fluidigm Flex Six IFC workflow (1006308) was employed for high-throughput amplification and target detection. Samples and assays were loaded into the IFC, followed by thermal cycling on the HX controller using the GE FlexSix Fast v2 protocol. Post-run data analysis was performed with Fluidigm Real-Time PCR software (68000088). GAPDH and HPRT1 served as endogenous controls, and gene expression was calculated using Ct values. All assays were run in duplicate.
Primary colonic fibroblast derivation and culture
Primary colonic fibroblasts were isolated and expanded from endoscopic biopsies obtained from healthy donors using the Lamina Propria Dissociation Kit (Miltenyi Biotec 130-097-410). In brief, biopsies were cut into 1–2 mm fragments and washed in DMEM, high glucose, GlutaMAX Supplement, pyruvate (Gibco #39166-021) supplemented with 10 mM HEPES (Gibco #15690-056) and Penicillin/Streptomycin (Sigma P0781). To remove epithelial cells, tissue fragments were incubated 3 times in HBSS (Gibco 14025-050) containing Pen/Strep, 10 mM HEPES, 5 mM EDTA (Thermo Fisher 15575-038), and 2 mM dithiothreitol (Thermo Fisher P2325) at 37 °C for 5 min each, with vortexing between incubations. This process was repeated 3 times in total (cumulative incubation ~45 min) and the remaining tissue fragments were used for stromal cell isolation.
For each digestion, tissue fragments were incubated in enzyme mix prepared according to the Lamina Propria Dissociation Kit instructions at 37 °C for 15 min, repeated 3 times. Between incubations, the tissue was mechanically dissociated using a blunt-end needle to facilitate digestion. The resulting suspension was filtered through a 70-µm cell strainer and centrifuged at 300g for 5 min at 4 °C. Pelleted cells were resuspended in DMEM supplemented with 10% FBS, Pen/Strep, and 1× ITS (Gibco 41400-045), and seeded into 6-well plates. Cultures were maintained at 37 °C with 5% CO2, and medium was changed every 2–3 days. Cells were passaged using TrypLE Express upon reaching confluency.
Lentivirus production
Lentiviral expression vectors for TWIST1 and OSR2 were obtained from F. Zhang via Addgene (1428908 and 144039, respectively)54. The empty vector was generated by excising the insert from the lentiviral vector using NheI-HF (NEB R3131S) and SpeI-HF (NEB R3133S) restriction enzymes, followed by gel purification and self-ligation of the linearized vector.
HEK293 cells were obtained directly from ATCC (CRL-1573), authenticated by the supplier, and tested mycoplasma-free. HEK293 cells were maintained in DMEM supplemented with 10% FBS and penicillin/streptomycin One day prior to transfection, 2 × 106 cells were seeded into T25 flasks. Transfection was performed the following day at 90–99% confluency. For each flask, 3.4 µg of the expression plasmid, 2.6 µg of psPAX2 (Addgene 12260), and 1.7 µg of pMD2.G (Addgene 12259) were transfected using 8.75 µl of Lipofectamine 3000 (Thermo Fisher L3000015), 7.5 µl of P3000 reagent, and 1.25 ml of Opti-MEM (Gibco 31985070). Culture medium was replaced 5 h post-transfection. Viral supernatant was collected 48 h later, filtered through a 0.45-µm PVDF membrane, aliquoted, and stored at –80 °C.
Lentiviral transduction
For transduction, primary fibroblasts were cultured to confluency in T25 flasks and incubated with an appropriate volume of lentiviral supernatant. After 48 h, the medium was replaced with fresh medium containing 1 µg ml−1 Puromycin (Thermo Fisher 15490717), which was refreshed every other day. Cells were passaged after seven days of selection. Lentiviral titres were estimated by transducing cells with three different volumes of virus and assessing cell viability following three days of complete Puromycin selection.
Immunofluorescence staining and quantification
Primary human colonic fibroblasts were cultured in 8-well µ-Slides (ibidi 80826). Following removal of the culture medium, cells were washed with PBS and fixed with 4% paraformaldehyde (Thermo Fisher J19943K2) for 30 min at 4 °C. Fixed cells were washed twice with PBS and permeabilized with 0.2% Triton X-100 (Sigma T8787) in PBS for 30 min at room temperature. After washing, cells were blocked in PBS containing 1% BSA (Cell Signaling Technology 9998) and 0.3% Triton X-100 for 30 min, then incubated at room temperature for 2 h with primary antibodies diluted in the same blocking solution: anti-Collagen VII (Abcam ab198899, 1:100) and anti-α-SMA (Dako M0851, 1:400). Cells were then washed with PBS 3 times and incubated for 1 h at room temperature in the dark with fluorophore-conjugated secondary antibodies, DAPI (Thermo Fisher D1306), and Alexa Fluor 647 Phalloidin (Thermo Fisher A22287). Cells were washed with PBS and mounted with Surgipath Micromount Mounting Medium (Leica 3801731). Samples were imaged with a Zeiss LSM980 confocal microscope and images were analysed using FIJI.
For quantification, approximately five images were acquired per condition to ensure sufficient cell numbers (>150 cells per condition). Maximum-intensity z-projections were generated for each fluorescence channel (Col7A, α-SMA and phalloidin), followed by thresholding and measurement of total signal intensity across the entire field of view. DAPI staining was used to identify nuclei and count cells using a standard FIJI workflow (contrast enhancement, Gaussian blur, thresholding, watershed segmentation, and particle analysis). Signal intensity for each marker was normalized to cell number, and for each patient, the normalized signal intensity in OSR2- or TWIST1-overexpressing conditions was calculated relative to the corresponding control.
Quantitative PCR
RNA was extracted from cultured primary fibroblasts using the RNeasy Plus Micro Kit (QIAGEN, 73034) and reverse-transcribed into cDNA using the High-Capacity cDNA Reverse Transcription Kit (Thermo Fisher 4374966). Quantitative PCR was performed using TaqMan assays (Thermo Fisher) according to the manufacturer’s instructions.
Bulk RNA-seq
RNA was extracted from cultured primary fibroblasts using the RNeasy Plus Micro Kit (QIAGEN 73034) and submitted to Novogene for bulk RNA sequencing. Library preparation and quality control were performed by Novogene, targeting 50 million 150-bp paired-end reads per sample.
Multiplexed immunofluorescence analysis
Section rehydration, antigen retrieval and blocking FFPE tissue sections were baked overnight at 60 °C and then rehydrated using the automatic slide processor Gemini Autostainer (Epredia/Shandon Diagnostics) by serial incubation for 5 min twice in each 100% xylene, 100% ethanol, 95% ethanol, 70% ethanol, 50% ethanol, and finally PBS. Subsequently, tissue slides were incubated in 0.3% Triton X-100 for 10 min and washed twice with PBS for 5 min. Antigens were retrieved using the NxGen decloaking chamber (Biocare Medical), wherein slides were heated to 110 °C for 4 min in pH 6 citrate buffer (S1699; Agilent or H-3300-250; Vector Laboratories) at 6.1 PSI then cooled to 100 °C over 16 min, rinsed in deionized water and transferred into a pH 9.0 Tris-EDTA buffer (12.1 g Tris base, 3.7 g EDTA, 1 ldistilled water) cooling from 100 °C to 85 °C over 20 min. Slides were then allowed to cool to room temperature and washed twice in PBS. Tissue sections were blocked overnight at 4 °C with 3% bovine serum albumin (A7906; Merck Life Science) and 10% donkey serum (C06SB; Bio-Rad) in PBS. Afterwards slides were washed twice with PBS, submerged in 1 mg l−1 of DAPI (D357; Thermo Fisher Scientific) in PBS for 15 min at room temperature and subsequently again washed twice in PBS. Mounting medium combining 50% glycerol (G5516; Merck Life Science) in PBS was then added to the sections and the slides covered with no. 1 glass coverslips (Leica Microsystems).
Iterative imaging, staining and dye deactivation slides were imaged using the Cell DIVE system (GE Research) through iterative rounds of staining with fluorescently labelled antibodies and deactivation of the fluorophores. Initially, entire slides were scanned to select the ROIs, then the selected regions were scanned to acquire the tissue-inherent autofluorescence signal. Subsequently, coverslips were removed by incubation in ice cold PBS and the sections incubated for 2 h at room temperature with three primary antibodies diluted in 3% bovine serum albumin and 10% donkey serum in PBS underneath a polycarbonate coverslip (Grace Bio-Labs). In the initial round primary antibodies were detected using fluorophore-labelled cross-adsorbed secondary antibodies raised in donkey (Thermo Fisher Scientific)—in all subsequent rounds primary antibodies directly conjugated to fluorophores were used. Fluorophores with excitation maxima around 488 nm, 555 nm and 647 nm were used, dependent on the antibody either Alexa Fluor 488, Alexa Fluor 555 and Alexa Fluor 647 (all Thermo Fisher Scientific) or CF 488 A, CF 555 and CF 647 (all Biotium). The stained slides were then washed twice in PBS, stained with DAPI and re-mounted as above. The selected ROIs per slide were then again scanned on the CellDIVE, and the autofluorescence signal previously acquired was removed from the signal of the staining round. The fluorophores were deactivated by de-coverslipping slides, washing them twice in PBS and incubating them at room temperature 3 times for 15 min in 0.1 M NaHCO3 at pH 11.2 (S6297; Merck Life Science) with 3% H2O2 (216763; Merck Life Science) with intermittent incubation in PBS for 5 min. After dye inactivation, slides were again incubated with DAPI, mounted and scanned as above to acquire the background anew. Subsequently, slides were stained with primary antibodies conjugated to fluorophores and imaged as in the first round. This background-staining-deactivation cycle was repeated until all markers had been acquired.
Primary antibodies were either purchased conjugated to fluorophores as indicated or conjugated with either the Mix-n-Stain Antibody CF488A/555/647 dye Labelling Kits (92233, 92234 and 92238; Biotium), the Alexa Fluor 488/555/647 Antibody Labelling Kits (A88062, A88065 and A88068; Thermo Fisher Scientific) or with a combination of the oYo-Link Thiol Antibody Reagent (AT3001; AlphaThera) and CF 488A maleimide, CF 555 maleimide and CF 647 maleimide (all Biotium). For the Mix-n-Stain and the Alexa Fluor kit conjugation followed manufacturer instructions, for the AlphaThera oYo-Link thiol kit, per 1 μg of antibody 33 pmol of oYo-Link thiol reagent and 155 pmol of maleimide-fluorophores were used. The oYo-Link thiol reagent was mixed with the maleimide-fluorophores in PBS and incubated at 37 °C on a heat block (Thermomixer comfort; Eppendorf) shaking for 2 h at 600 RPM in the dark. Subsequently, the resulting fluorophore-linked oYo reagent was mixed with the antibody and crosslinked to it for 2 h using 365 nm UV light at 4 °C (LED photo-cross-linking device AT8001-D; AlphaThera). Finally, one part of the conjugated antibody mix was diluted with two parts of PBS-based antibody stabilizer (131 050; Candor Bioscience).
Computational analysis
Raw sequencing data processing
All raw sequencing data were converted from BCL to FASTQ format using Illumina bcl2fastq (v.2.20.0.422) software. Up to one mismatch was allowed for each sample index barcode. Raw sequencing reads were then quality checked using FastQC software (v.0.11.9)55.
For each sequenced scRNA-seq pool, 10X Genomics Cellranger software (v.7.1.0) was used to process, align and summarize unique molecular identifier (UMI) and cell barcode counts into raw counts matrices. Human hg38 (refdata-gex-GRCh38-2020-A) reference genome was used for all alignments and gene annotation. Paired CITE-seq and hashing antibody panel data were processed together with gene expression libraries using the Cellranger feature barcoding workflow. Antibody UMI counts were summarized using custom made joint feature barcoding sequencing reference, pooling together TotalSeq antibody sequences of hashing and protein expression antibodies. Feature barcoding sequences are deposited at the Gene Expression Ominibus (GEO).
Visium FFPE ST library FASTQ files were similarly processed using 10X Genomics Spaceranger (v.2.1.0) software. Paired H&E images and Slide IDs used as input are deposited on GEO. Human hg38 reference genome was used, together with Visium Human Transcriptome Probe Set v.v1.0_GRCh38-2020-A.
scRNA-seq data analysis
Raw UMI counts matrices from gene expression and feature barcoding libraries were imported into R for further processing. Cell calling on raw gene expression matrices of all 10X whitelist barcodes was carried out using ‘emptyDrops’ function from DropletUtils R package (v.1.24)56. Raw counts matrices were further corrected for Illumina index swapping with swappedDrops57.
Cell QC metrics were calculated using scCustomize R package (v.2.1.2) and cell barcodes with low total UMI counts, low complexity and high mitochondrial RNA gene counts were filtered out from further analysis. The thresholds were determined per cell lineage basis from thresholding the empirical distributions after an initial clustering solution, as for instance epithelial cells typically contain higher percentage of mitochondrial UMIs than other cell types.
Following cell QC for each pool, cells were demultiplexed using hashing antibody counts matrices as follows. Non-hashing antibodies and tags absent from any given reaction were filtered out from the counts matrix. Counts matrices were normalized using centred log ratio transformation. Counts were clustered using clara k-mediods (k = number of samples + 1). The 99th percentile of a negative binomial distribution fit was defined as a positive threshold for each antibody tag. Doublets were identified as cells positive for multiple tags and filtered out from further analysis. Given sufficiently large number of samples in each pool, the majority of doublets will be heterotypic due to random mixing and therefore this approach enables robust, experimental label driven doublet removal. Cells which could not be confidently assigned to a sample label were also filtered out from further analysis. Demultiplexed cells were visualized as t-SNE embeddings and sex-specific gene expression (for example, XIST) was further examined to ensure the cells were segregated correctly as expected with sample-of-origin assignments and patient meta data.
All samples were merged using R package Seurat (v.5.1.0)58. Gene expression counts were normalized. Highly variable genes were detected separately within each sample pool in order to deemphasize between-batches variability and the union was used as input for dimensionality reduction by PCA. Principal components were then further batch-corrected using Harmony (v.1.2.1)59 algorithm for sample integration, and harmonized components were used as input for Louvain clustering and dimensionality reduction using uniform manifold approximation and projection (UMAP). Optimal clustering resolution was assessed using R package clustree (v.0.5.1)60 to assess cluster stability. Cell clusters were annotated using a combination of known marker genes and cross-classification with previously published scRNA-seq reference atlas datasets12,13,14,15,25,26,36,38,61 using Seurat label-transfer workflow.
scRNA-seq differential expression analysis
Differential gene expression analyses, including marker gene identification, were conducted using negative binomial generalized linear models for both ST and scRNA-seq data. To account for potential confounders, such as variability in gene detection rates and batch, donor or slide effects, these factors were incorporated as covariates into the model. Significance was assessed using Benjamini–Hochberg correction, considering genes with a false discovery rate (FDR) below 5% as significantly differentially expressed. Additionally, the MAST package62 was used to detect genes that exhibit more binary on/off expression patterns by modelling both continuous and discrete components.
For pseudobulk analyses of scRNA-seq data, raw UMI counts were aggregated to pseudobulk samples as previously described. The DESeq2 package63 was then used to normalize these counts and perform differential gene expression testing, again blocking for confounding variables. The Wald test was applied to derive P values, which were adjusted using the Benjamini–Hochberg correction.
Pseudobulk PCA analysis
Pseudobulk PCA analysis was carried out by aggregating UMI raw counts matrices of each individual sample following cell dehashing by summing counts for each gene. Epithelial, immune and stromal cells were considered separately. Sample pseudobulk counts were normalized using library size factor normalization, as implemented in DESeq2 R package (v.1.44.0)63 for bulk RNA-seq analysis. Counts were further transformed using variance stabilizing transformation (vst) and the top 1,000 most variable genes were selected to compute PCA. The first two principal components were visualized to assess the overall sample-level variability within each cohort. The above approach was also applied to assess sample-sample variability within the ST cohorts.
Differential abundance analysis
Cell-type or spatial region differential abundance analyses were carried out using three different approaches. R package miloR (v.2.0.0)64 was used to carry out graph-based differential abundance analyses based on data UMAP embeddings, an analysis which is independent of cell cluster assignments and can highlight differences between populations which may exist on a phenotypic continuum rather than falling into discrete clusters. In each case, integrated, batch-corrected harmony components were used for nearest-neighbour graph construction (k = 10). The R package sccomp65 (v.1.9.5) was used to further test cluster-level abundance differences. Cell proportions per sample were also compared using a two-sided Wilcoxon rank test.
Transcription factor regulon analysis
The pySCENIC pipeline (v.0.12.2b0)66 was utilized to identify active transcription factor modules in scRNA-seq datasets. Initially, the normalized single-cell gene expression matrix was filtered to exclude genes expressed in fewer than 20 cells. The RcisTarget database, which provides transcription factor motif scores for gene promoters and transcription start sites based on the hg38 human reference genome, was downloaded from the following resource. The expression matrix was then further filtered to include only genes present in the RcisTarget database.
Next, a gene–gene correlation matrix was constructed for co-expression module detection using the GENIE3 algorithm67. SCENIC66 was employed to analyse transcription factor networks, identifying co-expression modules enriched for target genes of candidate transcription factors from the RcisTarget database. Subsequently, the AUCell package was used to calculate a score for each transcription factor module in individual cells.
To pinpoint condition- or cluster-specific transcription factor modules, generalized linear models were applied to test for the dependence of transcription factor AUC values on conditions or clusters, incorporating batch effects and gene detection rates as covariates in the model to account for the significant influence of detection rate on AUC values. The resulting P values were adjusted for multiple comparisons using the Benjamini–Hochberg correction method.
Cell cycle scoring
Cell cycle phases were predicted using Seurat’s CellCycleScoring function, using human 2019 cell cycle gene reference.
Xenium ST cell segmentation
Cell segmentation for Xenium ST was conducted by first re-segmenting cells using a transcript-density-based approach implemented in Baysor algorithm68. Nuclei-based segmentation from XeniumRanger outputs was integrated as a prior with a weighting factor of 0.7. –n-clusters parameter was set to 12, determined from the initial clustering solution of nuclei-segmented data. Following segmentation, a new cell-by-gene matrix was constructed by including only transcripts with a transcript assignment probability above 0.9. Transcripts with a confidence score below 0.9 were subsequently excluded from all downstream analyses. These analyses were similarly repeated for Xenium Prime assay, except cell boundary image-based cell segmentation was used as a prior instead.
Cell segmentation quality was assessed both visually and via transcriptome metrics. We calculated silhouette scores, as well as mutually exclusive co-expression ratio (MERC)69 metrics to evaluate the amount of transcript wrongly assigned from adjacent cells. For MERC score, we assessed both universal cell-type markers, as well as scores focused on problem areas in intestinal tissue specifically. These included T cells and B cells, which form dense aggregates in follicles and therefore are difficult to segment and epithelial cells and T cells, where intra-epithelial lymphocytes wedged in between epithelial cells are often segmented as part of the nearby epithelial cells. MERC gene sets used are provided as part of Mendeley Data.
Xenium ST clustering analysis
Cell by gene counts matrix following cell re-segmentation was converted into Seurat R object for further analysis. Very small cells, which likely represent segmentation artefacts, and cells with very few detected transcripts (<15 per cell) were filtered out from further analysis. Cells with higher than expected (>99th percentile) negative probe signal were also filtered out from further analysis. As the remaining cells formed a prohibitively large dataset (>7M cells in total), undertook a sketch-based analysis workflow, as implemented in Seurat R package. In brief, for each tissue section, we sampled 5,000 cells using the LeverageScore method to select most informative cells. Data from the sketched subset cells was then dimensionality-reduced using PCA, followed by slide integration using Harmony, as described before. Nearest-neighbour graph was constructed for downstream clustering analysis and clusters were visualized using UMAP embeddings. The remaining cells from the full dataset were then projected onto the sketched data and onto the sketched UMAP model representation using ProjectIntegration and ProjectData functions. All subsequent analyses and statistics were computed on integrated dataset from all slides.
Low resolution clusters were annotated as broad cell-type lineages based on marker gene expression and label transfer of scRNA-seq reference data onto the sketched dataset. Each broad lineage identified (such as epithelium, fibroblasts or T cells) was then subset and further sub-clustered to identify cell subsets, which were similarly annotated as before. In the case where clusters were identified by our Xenium panel markers that did not have a one-to-one relationship to a scRNA-seq reference cluster and/or our panel did not contain subtype-defining markers, we annotated clusters based on their cell lineage combined with most specific gene marker expression (for example, SPP1+ macrophages, NOS2+ epithelium).
Xenium ST cell niche detection
Cell niche detection was performed using an unsupervised approach implemented as part of the spamsc R package (https://github.com/ChloeHJ/spamsc). Initially, for each field of view within the ST dataset, the spatial coordinates of cell centroids were extracted to determine the nearest neighbours for each cell, with the number of neighbours set to 30. Utilizing these coordinates, two separate neighbour matrices were constructed for each cell: (1) expression counts of each cell’s nearest neighbours were aggregated, generating a cell-by-gene matrix that encapsulates the local transcriptomic environment; and (2) cell types, as defined by high-resolution cell subclusters from previous analyses, in each cell neighbourhood were counted and aggregated into a cell-by-cell-type matrix that encapsulates the local cell-type neighbours for each cell. This process was iteratively applied across all slides and the matrices were merged. To identify tissue niches, unsupervised clustering of the neighbourhood matrices was used by first carrying out PCA (for gene expression matrices only), integrating the principal components using Harmony across different slides and detecting niche clusters using Louvain clustering. Niches were labelled according to their location, tissue architecture features or cell-type enrichment.
Xenium ST cell co-localization analysis
For cell–cell co-localization analysis in Xenium ST dataset, the following approach was taken. First, using spatial coordinates of each cell centroid, nearest 30 spatial neighbours for each cell were determined and cell-type counts in each spatial cellular neighbourhood were aggregated. We then constructed a cell-type correlation matrix between each cell and its neighbours, representing the frequency of co-occurrence between different cell types. A separate co-occurrence matrix was constructed for each different experimental group. Cell-type co-occurrence networks were visualized as graphs, where edges of low co-occurrence (<0.15), as well as negative co-occurrence (that is, cells always localize away from each other—for example, muscle cells and epithelium) were filtered out. Graphs were laid out using force-directed layout and visualized using ggpraph R package (v.2.2.1).
Xenium ST custom panel design
A custom 480-gene panel was designed for Xenium ST to delineate over 100 fine-grained cell types within the gut, as defined by previously published single-cell reference datasets and scRNA-seq data generated as part of this study. Additional genes from literature were included to cover cell types that are typically difficult to capture using droplet-based scRNA-seq approaches—for example, neurons and neutrophils, as well as squamous epithelium markers, a cell type that is typically not present in the intestine. Genes were expert-selected to ensure comprehensive representation across stromal, immune, and epithelial lineages while excluding highly abundant transcripts to prevent optical crowding. For each identified cell type, 5 to 10 redundant genes were incorporated to ensure cell typing robustness. The panel’s performance was iteratively evaluated by assessing its ability to reconstruct scRNA-seq clustering patterns in the reference datasets. Probe balance was further optimized using 10X online Xenium Panel Designer tool. Xenium Prime 100 gene custom add-on panel was designed using the same procedure, focusing on key missing intestinal marker genes and disease-specific markers. All custom panel designs have been provided as part of Mendeley Data.
scRNA-seq meta-analysis
For scRNA-seq meta-analysis, data were downloaded from GEO, ENA or Zenodo repositories (see Data availability). Individual datasets were initially clustered to replicate analyses from their respective studies and fibroblast clusters from each study were subset for further downstream analyses. Fibroblast cells from public datasets and this study were then merged and integrated using harmony, correcting for both between dataset differences and known within-dataset batch variables.
To identify transcriptional programmes across transcriptomic data, cNMF70 was applied to gene expression matrices. Prior to factorization, harmony59 was used to correct for donor-specific batch effects. cNMF was run to extract 30 factors, each representing a gene expression programme recurrent across samples. Resulting gene loadings and usage scores were used for downstream interpretation and spatial mapping of inferred programmes. To visualize the relative contribution of cNMF-derived gene expression programmes (GEPs) across cell states, we calculated per-cluster usage proportions. First, harmony-corrected cNMF usage scores (that is, per cell factor loadings) were integrated into the Seurat metadata. For each annotated cell cluster, average usage scores were computed per factor. GEPs with low average usage (<0.05) were excluded to reduce noise. Remaining average scores were normalized within each cluster to sum to one, yielding relative proportions per GEP. These proportions were visualized as pie charts, with all non-disease cluster associated GEPS collapsed into an ‘other’ category for ease of visualization. All visualizations were generated using ggplot2 with coord_polar() for circular layout.
cNMF factors derived from meta-analysis were then applied to Xenium in situ 5100-plex dataset, calculating activity scores for each cell in the spatial dataset. Furthermore, cells from a diabetic wound-healing dataset37 were similarly scored. Separately, cNMF factors were similarly calculated for fibroblast cells from Xenium in situ 5100-plex dataset, and reciprocally applied to the scRNA-seq meta-analysis to cross-map gene expression programmes identified with these different approaches.
Xenium spatial autocorrelation analysis
To identify gene pairs with spatially correlated expression patterns, we used two approaches. Firstly, local spatial neighbourhood gene expression matrices were computed as described in ‘Niche analysis’ and gene–gene Pearson’s correlation was calculated to detect co-localization gene expression patterns regardless of cell type. Additionally, for each gene pair, we also calculated a bivariate Moran’s I statistic, which quantifies how the spatial distribution of each gene related to that of another. Statistical significance was assessed by permutation testing, in which the spatial positions were repeatedly randomized to generate a null distribution of Moran’s I values. The resulting P values were adjusted for multiple testing, and significant positive or negative bivariate Moran’s I values indicated spatially co-expressed or contrasting gene pairs, respectively.
scRNA-seq spatial projection
Spatial projection of fibroblast scRNA-seq data was carried out using the spamsc R package workflow (https://github.com/ChloeHJ/spamsc). In brief, to project single-cell fibroblast data onto a Xenium spatial dataset, we first identified genes shared between the scRNA-seq fibroblast dataset and the Xenium data, removing any features not expressed in either dataset. After merging these filtered datasets into a single Seurat object, we normalized and scaled them, then applied Harmony to batch correct for differences between the single-cell and spatial modalities. We extracted a reduced dimensional embedding and used a k nearest-neighbour approach to map fibroblast single cells onto the Xenium slide. For each fibroblast cell in the single-cell dataset, we searched for neighbouring Xenium spatial cells in the integrated reduced space; if spatial fibroblast cells were among the top neighbours, the spatial coordinates of the closest such neighbour were assigned to that cell. The accuracy of this mapping was assessed by comparing gene expression patterns between projected fibroblasts and the local spatial environment, using correlation analyses and evaluating enrichment of known fibroblast markers. These methods were applied using SPAMsc functions MergeDatasets and RunProjection.
Xenium cell adjacency-dependent gene expression analysis
To assess how a given cell type’s gene expression patterns differ depending on the identity of its neighbouring cells, we implemented a procedure that aggregates gene expression counts based on spatial adjacency. First, we extracted cell centroid x,y coordinates for each tissue section. For each section, a nearest-neighbour graph was constructed using Euclidean distances between cell centroids, and for each cell, the k nearest neighbours were identified. We then focused on a particular cell type or group of cell types—for example, only epithelial cells, or FAS cells. Cells not belonging to the target identity group were excluded by zeroing out their counts. The expression values of the selected gene set (from a specified spatial assay) for the k nearest neighbours of each cell were then summed, yielding a ‘niche expression matrix’ that represents the aggregated gene expression within local neighbourhoods populated by the specified cell types. By applying this process across multiple images and combining the results, we obtained a matrix that captures how gene expression in a cell of interest changes when it is surrounded by different cell populations. The aggregated niche expression matrix was then filtered to remove cells which do not have any cells of the selected cell type in their local neighbourhood. Gene expression counts of the remaining cells were then normalized and subjected to downstream differential expression analyses to quantify gene expression differences conditioned on spatial adjacency to particular cell types.
Distance gradient-based analyses
For fistula edge gradient-based analyses, fistula were selected where the tract was clearly delineated, with samples excluded where the edge of the tract was not fully captured. Fistula tract edges were then manually annotated using interactive XeniumExplorer software and coordinates of the polygonal annotations were exported as csv files. These were imported and converted into sf polygon objects. Then for each tissue section, centroid coordinates of tissue-detected cells were extracted, converted into spatial features, and distances to the nearest polygon edge were computed using st_distance() from the sf package. These distances were added to the Seurat metadata and visualized using ImageFeaturePlot().
For downstream analysis, cells were grouped by cell type and donor (slide ID), and kernel density estimates of cell abundance as a function of distance were computed independently per donor using the density() function in R. To summarize across donors, mean densities were calculated at each distance point, and standard errors were estimated as the standard deviation of densities across donors divided by the square root of the number of samples per group. These were visualized as mean ± standard error using geom_line() and geom_ribbon(). To identify genes with expression patterns that vary as a function of distance from fistula edges, we modelled gene expression counts using negative binomial generalized additive models with distance as a smooth spline term. For each gene, significance of the distance effect was assessed by comparing the full spline model to a null model using likelihood ratio test. These tests were carried out within lineage—that is, distance varying genes were identified in fibroblasts and macrophages separately by subsetting these cells individually. Genes with significant spatial associations were visualized using heat maps, plotting loess-smoothed gene expression over distance.
Epithelial-edge distance-based analyses were carried out as above, except starting from manually annotated epithelial edges and excluding all non-epithelial cells from the analysis.
Picrosirius red image analysis
Individual images were first processed in Qupath software (v.0.5.0). In brief, tissue-covered areas of each slide slides were first selected as ROIs, with areas containing imaging or tissue artefacts filtered out. A groovy script was then used to systematically partition and export image tiles of defined size (512 × 512 pixels) with partial overlap (64 pixels), with x,y coordinates for each image tile exported along the image for downstream processing.
All subsequent image analyses were then performed in Python (v.3.9). All analyses were conducted using skimage for image processing71,72, numpy and pandas for data handling, matplotlib and seaborn for visualization, and scipy and scikit-learn for statistical and machine learning operations. For each image, the red colour channel was extracted, pixel intensities were normalized to 0–1 range and uniform global thresholds of 0.15 and 0.18 for all images were set for red channel binarization, after initial exploratory analyses with Otsu’s thresholding to guide parameter settings. The red channel image was then binarized at the two thresholds for parallel analyses, with the higher threshold excluding less well stained, sparser collagen bundles. Following binarization, we undertook morphological pre-processing of each image tile with erosion and opening using skimage.morphology operations to remove small specs of noise and smooth the collagen structures. Following pre-processing, each binarized tile was skeletonized using skimage.morphology.skeletonize, producing one-pixel-wide paths representing the collagen bundle backbone structures. Small, spurious branches were pruned. The resulting skeleton was used to compute several morphological and structural metrics as follows.
-
1.
Length and thickness: connected components in the skeleton were first identified using skimage.measure.label and their perimeter approximations were used as a proxy for collagen length distribution. This yielded total length, mean length and selected quantiles (5th, 50th and 95th) as features for each tile. Thickness was estimated by calculating Euclidean distance transform with scipy.ndimage.distance_transform_edt on the binary mask and extracting the distance values along the collagen skeleton. Total, mean and quantile values for each image patch were calculated as features.
-
2.
Connectivity and morphological descriptors: the number of connected components within the skeleton was used to calculate collagen network connectivity in each image tile. We also calculated perimeter-area ratio, eccentricity (shape elongation) and average elongation (major/minor axis ratio of fitted ellipses) using skimage.measure.regionprops.
-
3.
Fractal and textural measures: fractal dimension was estimated using a box-counting method across multiple scales. Lacunarity, a measure of textural heterogeneity, or ‘gappyness’, was computed at four different box sizes for each image patch. Additionally, Shannon entropy, which can provide a measure of image complexity, was calculated using skimage.measure.shannon_entropy.
-
4.
Anisotropy and orientation: anisotropy was calculated from the image patch skeletons by calculating the variance in the orientations of the skeleton segments. Anisotropy measures the overall directional coherence of the network; that is, how parallel each segment is to each other.
-
5.
Mean free path and collagen fraction: a mean free path along both collagen regions and regions between collagen areas was estimated by computing pairwise distances between feature points. To compute the total collagen fraction, the proportion of thresholded pixels at both thresholds used was calculated for each image path.
All metrics were mapped to x,y spatial coordinates of the image patches for further spatial analysis and visualization.
For image clustering analysis, a variational autoencoder was trained on exported image tiles, integrated with a pretrained VGG16 encoder as a feature extraction backbone. In brief, a custom image tile generator class was implemented to load and preprocess image tiles to resize, normalize intensities and apply gamma correction. The generator also computed Shannon entropy for each tile to determine sample weights to reduce the weights of any texturally uninformative tiles (for example, edge tiles) during training. For VAE architecture, the encoder consistent of an initially frozen VGG16 network pretrained on ImageNet to extract robust feature representations, followed by a fully connected layer and two dense layers producing the latent parameters with a latent dimensionality of 64. A sampling layer generated latent vectors from these parameters. The decoder employed a series of transposed convolutional and up sampling layers to reconstruct the original image resolution from the latent space. The entire model was compiled with a mean-squared-error loss and trained until convergence on total of 301,941 image tiles. 64 dimension latent vector was then extracted for each image tile and used as input for Louvain clustering (resolution – 0.3). Clusters were intersected with morphological and structural tile metrics calculated as described above, and annotated on the basis of average collagen density per cluster, with highly fibrotic tissue regions annotated as ‘very high’ clusters and tissue regions with lower levels of collagen deposition (such as mucosa) as ‘low’ or ‘mid’ regions. Comparative analysis of relative cluster abundance between fistulae and control samples were carried out using sccomp R package, as described above.
Fibroblast lineage prediction analysis
pySCENIC was used to reconstruct transcription factor regulons for fibroblast cell subsets specifically and cell-by-AUC matrices for each transcription factor were imported into R for further processing. Transcription factor activity scores were scaled and used as input for dimensionality reduction by PCA. Fibroblasts cells from healthy control samples were grouped into broad fibroblast subtypes -Stromal 1–4 subsets and served as a reference dataset, while all other cells, including fistula-associated fibroblasts, were then classified with respect to the reference using Seurat’s label-transfer workflow (FindTransferAnchors and TransferData). For each query cell, a prediction probability was calculated on the basis of transcription factor regulon activity score embeddings, scoring the cell as most similar to each of the four main fibroblast subtypes, Stromal 1–4.
Pathway enrichment analysis
Gene Ontology (GO) and Reactome pathway enrichment analyses for both scRNA-seq and ST data were conducted using the clusterProfiler R package73 (v.4.12.6). Gene identifiers were mapped with the org.Hs.eg.db (v.3.19.1) annotation DBI package (v.1.66.0). Individual cluster marker sets and differentially expressed genes were analysed for overrepresentation, using all expressed or detected genes as the background reference. Hypergeometric P values were adjusted for multiple testing through the Benjamini–Hochberg method. The enrichment results were visualized with clusterProfiler and ggplot2 packages. For pathway activity analysis at the single-cell or spot level, pathway scores were derived as gene module scores using the AddModuleScore function in the Seurat R package.
Visium ST clustering analysis
Raw UMI count matrices output by 10X Spaceranger pipeline, corresponding images, spot coordinates, and scale factors were first imported into R for further processing. The spot matrix was filtered to retain only spots that overlapped with tissue. To assess background signal, we fitted a negative binomial distribution to the total UMI counts of spots located outside tissue sections. Using this distribution as a reference, we further removed any tissue-overlapping spots whose UMI counts were indistinguishable from the non-tissue background, as these likely represented areas of under-permeabilization. Most spots removed via this procedure were either section-specific anomalies or associated with tissue artefacts. Additional tissue irregularities (e.g., folds, tears) were annotated based on H&E images, and the corresponding spots were also excluded from the downstream analysis manually. Spots under coverslip air bubbles, however, were retained, as these defects affected only the H&E image and not the underlying tissue integrity.
Next, raw UMI counts were normalized using SCTransform, a regularized negative binomial regression-based method74, to account for variability in total RNA content across spots. PCA was used for dimensionality reduction, using the union of spatially variable genes across all slides detected using Moran’s I spatial autocorrelation method. For integrative analysis, spots from individual slides were combined using the Harmony59 algorithm to correct for slide-specific effects. Subsequent clustering analysis was performed using the Louvain algorithm at a resolution of 0.5, and results were visualized using UMAP. Clusters were visualized on H&E images with a spot size scaling factor of 1.6. The resulting integrated clusters were compared to those derived from individual slides to confirm that no biologically relevant heterogeneity was lost after batch correction. In addition, we ensured that corresponding tissue regions across slides converged into the same integrated clusters, confirming successful alignment of equivalent anatomical areas.
Visium ST cell-type deconvolution
Visium ST cell-type deconvolution was carried out using R package CARD75. Transcriptomic signatures and histopathological features detected in our Visium ST analysis indicated the presence of both squamous epithelium and neutrophils in our disease tissue sections, which were not captured by our scRNA-seq analysis. To account for this missing data and ensure a more accurate, robust cell-type deconvolution result, a combined scRNA-seq single-cell reference dataset was constructed as follows. scRNA-seq from the intestine generated in this study was merged with neutrophil cells identified in the GEO GSE163668 dataset76, as well as squamous epithelium cells identified in the GEO GSE201153 dataset77, the data were normalized, merged and re-embedded together using only common genes between gene matrices of all three reference datasets and Visium ST data.
Visium ST cell–cell signalling analysis
To identify region-specific, spatially co-localized cellular signalling events, we undertook analysis as previously described in36. We first obtained receptor–ligand databases from78 and79. Each ST spot was scored for receptor–ligand co-expression by incorporating weighted expression from neighbouring spots. We calculated pairwise Euclidean distances for all ST spots and assigned distance-based weights, setting weights to zero for spots beyond two neighbours and normalizing weights for closer spots. For each spot, we computed a distance-weighted receptor–ligand product score, scaled by the number of contributing spots to account for tissue edges.
To determine significance, we performed 100 random permutations of spot locations to generate an empirical background distribution and calculated P values for each receptor–ligand pair based on this distribution. After applying Benjamini–Hochberg correction for multiple testing, spots with FDR < 5% were considered positive for receptor–ligand cross-talk. We then prioritized region-specific interactions using generalized linear modelling, assessing the dependence of receptor–ligand scores on spatial clusters while controlling for gene detection rates. Condition-specific interactions were similarly modelled, and all significant receptor–ligand pairs were further subject to multiple testing correction to ensure robust identification of region- and condition-specific cellular communication events.
Bulk RNA-seq analysis
Raw sequence reads were quality checked using FastQC software55. Cutadapt80 was used to trim poor-quality bases and Illumina universal adapter sequences from raw reads before alignment.
The human hg38 reference genome analysis set was obtained from the University of California Santa Cruz (UCSC) ftp site81. Full lentiviral construct sequences (Addgene plasmids #142908 and #142826) were included as additional reference contigs and indexed together with hg38 as a reference genome using STAR aligner82. Reads were then aligned to this custom reference.
Picard tools83 were used to mark duplicate sequences as an additional quality control step. Raw gene expression counts were summarized with featureCounts84 using a custom hg38 and plasmid GTF file containing joint annotations. The MultiQC tool was used to aggregate quality metrics. Sample quality metrics and raw read counts were imported into R for further processing. The DESeq2 R package63 was used to estimate library size factors, normalize counts and perform differential expression analyses. Benjamini–Hochberg multiple testing correction was used to compute FDR, and genes were considered significantly differentially expressed at <5% FDR. PCA was performed in R using the top 1,000 most variable genes, with normalized DESeq2 variance-stabilized transformation expression as input. Combat85 was used to correct donor-specific effects for heat map visualizations only.
Multiplexed immunofluorescence quantification
Multiplexed immunofluorescence images were analysed by selecting ROIs in each sample, avoiding areas with tissue artefacts or mismatched initial and final DAPI staining. In fistula samples, ROIs were annotated as granulation tissue, lesion, lesion edge or fibrotic zone, in order to be consistent with ST region definitions. Cell segmentation was performed in QuPath software (v.0.5.0) using DAPI for nuclear detection and expansion of cell boundaries by 5 µm to approximate full cell outlines. Mean marker intensity was then quantified per cell. To enable cross-sample comparison, intensities were clipped at the 95th percentile within each sample to exclude high-intensity outliers. Cells were binarized as positive or negative per marker based on sample-specific thresholds, and the proportion of marker-positive cells was calculated per region and compared across samples.
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
