Patient recruitment and sample acquisition

CTCL skin samples generated for the present study were donated with written consent and approval from the Newcastle and North Tyneside NHS Health Authority Joint Ethics Committee (protocol 08/H0906/95+5). Each patient donated two skin-punch biopsies, from a representative plaque or tumor. One biopsy was used for scRNA-seq and the other for bulk sequencing and IHC. All patients had mycosis fungoides, diagnosed based on correlation of clinical and histopathological features. The stage of CTCL at the time of the biopsy was based on clinical assessment performed by dermatology specialists at the Department of Dermatology and National Institute for Health and Care Research (NIHR) Newcastle Biomedical Research Centre. For fixed sample acquisition for both IHC and single-cell flex, FFPE blocks were accessed through the Newcastle CEPA Biobank (17/NE/0070). For additional IHC validation cohorts, samples were donated with consent from the local ethics committee at the Medical University of Vienna (ECS 1360/2018) and the Swedish Ethical Review Authority (2019-03467). CTCL diagnosis and staging as well as monitoring for disease progression were performed by specialists in dermatology and dermatohistopathology at the Department of Dermatology, Medical University of Vienna and the Department of Dermatology, Karolinska University Hospital, Stockholm. Detailed patient information including sex and age is provided in Supplementary Table 1.

Fresh skin sample processing

Skin biopsies in eight patients were immediately processed by removing the lower dermis and subcutis and separating the epidermis and dermis after dispase II digestion at a concentration of 2 U ml−1 for 2–3 h at 37 °C. The epidermis and dermis were processed separately in type IV collagenase at a concentration of 1.6 mg ml−1 overnight (37 °C, 5% CO2). Subsequently, single-cell suspensions were formed by vigorous pipetting and filtering (100-μm filter), counted and further processed via FACS.

FACS sorting and 10x Genomics Chromium loading

Cells from both the epidermis and the dermis were stained with an antibody panel containing CD45 (BUV395, clone HI30, BD Biosciences, 563792), CD8a (Alexa Fluor-700, clone HIT8a, BioLegend, 300920) and sorted using FACS into the following fractions: CD45−, CD45+CD8a+ and CD45+CD8a−. FlowJo was used to analyze FACS data. A target of 10,000 cells was used to calculate the loading volume for the 10x Chromium, taken from the manufacturer’s protocol (10x Genomics User Guide, CG000207). Each fraction sorted from the epidermis and dermis was loaded on to one channel of the 10x Chromium chip before running on the Chromium Controller using the Chromium Next GEM Single Cell 5′ v.1.1 kit.

Library preparation and sequencing

Gene expression libraries were generated from the resulting complementary DNA after clean-up following the 10x Genomics protocols. Enriched TCR cDNA was also generated from each CD45+ fraction (10x Genomics Chromium Single Cell V(D)J Enrichment Kit, Human T Cell) and subsequent libraries were made. All libraries were sequenced using an Illumina NovaSeq with the gene expression libraries sequenced to achieve a minimum of 50,000 reads per cell and the T libraries sequenced to achieve a minimum of 5,000 reads per cell.

Single-cell fixed RNA profiling

CTCL FFPE sample blocks were utilized for single-cell fixed RNA profiling using the single-cell gene expression flex kit from 10x Genomics. Sixteen samples from fourteen patients were sectioned using a microtome (Leica, RM2235) and, for each sample, four sections were cut at 25-μm thickness. Sample sections were deparaffinized and pestle dissociated according to the manufacturer’s guidelines (10x Genomics User Guide, CG000632) and stored short term at 4 °C. Four sample pools, each containing four samples, were generated by multiplexing and probe hybridized for 18 h, followed by gel beads-in-emulsion generation, barcoding and library preparation according to the 10x Genomics user guide (CG000527) and sequenced down one lane of a Novaseq 6000 S4 flow cell.

Whole-genome sequencing

DNA was extracted from frozen skin samples in RNA Later (Invitrogen) using the AllPrep Micro kit (Qiagen, 80284) following the manufacturer’s protocol. The DNA was quantified using a Qubit with the High Sensitivity DNA kit (Invitrogen, Q32851). Library preparation was carried out using NEBNext Ultra II DNA Library Prep Kit from Illumina. Libraries were uniquely dual indexed to mitigate tag hopping, quantified, equimolar pooled and sequenced on an Illumina NovaSeq 6000 platform (150-bp paired-end reads).

FFPE Visium CytAssist spatial transcriptomics

RNA quality and tissue morphology of CTCL FFPE sample blocks were assessed before FFPE Visium processing. Each sample block was sectioned using a microtome (Leica, RM2235) at 5-μm thickness on to a SuperFrost Plus microscope slide (VWR, 6310108), incubated for 3 h at 42 °C, dried overnight in a desiccator at room temperature and processed for FFPE Visium within 2 weeks of sectioning. Deparaffinization, hematoxylin and eosin (H&E) staining and decrosslinking steps were performed as per the manufacturer’s recommendations (10x Genomics Demonstrated Protocol, CG000520) and sections were imaged on a Hamamatsu Nanozoomer. Sections were then further processed with FFPE Visium CytAssist v.2 chemistry (6.5-mm) kit and dual-indexed libraries were prepared as per the 10x Genomics User Guide (CG000495). Four libraries were pooled at a time and sequenced down one lane of an Illumina Novaseq SP flow cell with the following run parameters: read 1: 28 cycles; i7 index: 10 cycles; i5 index: 10 cycles; read 2S: 50 cycles.

Fresh frozen Visium spatial transcriptomics

Healthy skin samples were embedded fresh into Optimal Cutting Temperature embedding medium and frozen using an isopentane-dry ice slurry. All tissues were sectioned at a thickness of 15 µm and a tissue optimization experiment using Visium Spatial Gene Expression Reagent Kits—Tissue Optimization (10x Genomics) was carried out to determine the optimum permeabilization time (14 min) using the manufacturer’s protocol. Spatial gene expression libraries were then generated using the Visium Spatial Gene Expression Reagent Kits (10x Genomics) and the manufacturer’s protocol and sequenced using an Illumina NovaSeq 6000 to achieve a minimum number of 50,000 read pairs per tissue-covered spot. H&E images needed for analysis were taken using a Zeiss AxioImager with an apotome microscope (Carl Zeiss Microscopy) and Brightfield imaging (Zeiss Axiocam 105 48-color camera module) at ×20 magnification. The ZEN blue edition v.3.1 (Carl Zeiss Microscopy) software was then used to acquire the images and adjust the z-plane and light balance, as well as stitching the image tiles to retrieve the overall H&E image file.

RareCyte 16-plex immunofluorescence staining

FFPE sample blocks were sectioned using a microtome (Leica RM2235) at 5-µm thickness. Slides were dried at 60 °C for 60 min to ensure adherence. After deparaffinization, sections were subjected to antigen retrieval for 5 min at 95 °C followed by 5 min at 107 °C using EZ-AR 2 Elegance AR buffer (pH 8.5; BioGenex, HK547-XAK). To remove autofluorescence, slides were bleached with AF Quench Buffer which consists of 4.5% H2O2/24 mM NaOH in phosphate-buffered saline (PBS) and quenched for 60 min with a strong white light exposure (using the HIGH setting) followed by further quenching for 30 min using an ultraviolet transilluminator (365-nm HIGH setting). Slides were rinsed with 1× PBS and incubated in 300 µl of Image-iT FX Signal Enhancer (Thermo Fisher Scientific, I36933) for 15 min. After rinsing with surfactant wash buffer (0.025% Triton X-100 in 1× PBS), the slides were stained with 300 µl of labeled primary antibodies in mouse/rabbit diluent (PBS antibody stabilizer (CANDOR Bioscience, 131500) containing 5% mouse serum and 5% rabbit serum) for 120 min. All antibodies (Supplementary Table 8) were pre-diluted according to the manufacturer’s recommendations. Slides were washed with a surfactant wash buffer and incubated with 300 µl of Hoechst 33342 (Thermo Fisher Scientific, H3570) in goat diluent (10% goat serum in PBS antibody stabilizer) for 30 min, after washing in 1× PBS. Finally, slides were coverslipped using ArgoFluor mounting medium (RareCyte, 241301000) and imaged on the following day using a RareCyte Orion microscope with a ×20 objective. Scans were performed using Imager and relevant acquisition settings were applied using the software Artemis.

scRNA-seq data processing

We processed newly generated and published raw scRNA-seq data using a consistent pipeline. In brief, raw sequencing data for published datasets (National Center for Biotechnology Information (NCBI) accession HRA000166 and Sequence Read Archive (SRA) accession PRJNA754592) were downloaded from the National Genomics Data Center (NGDC) and the SRA, respectively. For HRA000166, we downloaded raw data for both gene expression and V(D)J because they are publicly available. For accession PRJNA754592, only gene expression raw data were available, because the authors of the original study submitted only bam files from which V(D)J fastqs cannot be restored. As a result, V(D)J data from accession PRJNA754592 were not included. Gene expression data from droplet-based sequencing were processed using STARsolo64 (v.2.7.10a_alpha_220818) and the 2020-A reference genome. STARsolo outputs were then processed with cellbender65 (v.0.2.1) to remove ambient RNA contamination. For the flex samples, CellRanger (v.8.0.0) was used for cell calling with Chromium_Human_Transcriptome_Probe_Set_v1.0.1_GRCh38-2020-A as a reference. We initially had 16 samples from 14 patients. Samples with <400 cells were removed, so there were finally 11 flex samples from 10 patients. V(D)J data were processed using the 10x software package CellRanger vdj (v.7.2.0) and the ‘refdata-cellranger-vdj-GRCh38-alts-ensembl-7.1.0′ reference genome. Gene expression outputs from STARsolo were read using the read_10x_mtx function in Scanpy66 (v.1.8.1 and v.1.9.3). Data objects from different 10× lanes were then concatenated using the concatenate function in anndata (v.0.7.6). To detect and remove doublets, we applied Scrublet67 (v.0.2.3) to the data from each 10× lane to obtain per-cell scrublet scores and used a doublet exclusion threshold of median plus four median absolute deviations of the doublet score, as previously described22. Cells with >20% mitochondrial gene expression or expression of <400 detected genes, or total gene counts <1,000, were excluded from downstream analysis. Genes that were expressed in fewer than three cells were also removed.

Data clustering and integration

After preprocessing and quality control of the data, we created two major data objects: (1) CTCL plus skin cell atlas object, and (2) CTCL-only object through data integration. To create the first data object, we integrated the newly generated CTCL scRNA-seq data (from both fresh and FFPE samples) with CTCL data from accessions HRA000166 and PRJNA754592 and healthy skin, AD and psoriasis data from an existing skin cell atlas22. The integration was carried out using the scVI module in scvi-tools68 (v.0.20.3) as the following steps. The highly variable genes were selected by sc.pp.highly_variable_genes (n_top_genes=4000, batch_key=‘donor’, span=1). The model was built with (n_hidden=100, n_layers=1, n_latent=6, gene_likelihood=‘zinb’), using batch_key=’donor’ and categorical_covariate_keys= [‘technology’] (10x versus Flex). The integrated object was clustered by Leiden clustering at a resolution of 1. To broadly annotate cell types, gene expression counts were normalized and log(transformed). The expression of main cell markers was visualized on the main object and clusters were annotated accordingly. For the granular annotation, we used CellTypist to train a classifier model based on annotations from the skin cell atlas22 and then used the model to annotate the integrated object. However, as the model does not include B cells because they are not in the skin cell atlas data, the B cell cluster was manually annotated based on marker gene expression. To make the second data object, we integrated the newly generated CTCL scRNA-seq data (from both fresh and FFPE samples) with CTCL data from accessions HRA000166 and PRJNA754592. The top 4,000 highly variable genes among donors were selected using sc.pp.highly_variable_genes function in scanpy. Then scVI model was set up with donor as batch key and study and technology as technical covariates. The following parameters were used for the model: n_hidden=100, n_layers=1, n_latent=5, gene_likelihood=‘zinb’. As all CTCL samples are already part of the first object (with skin atlas samples), annotations were transferred from that object to the CTCL object.

Differential abundance analysis using Milo

To reveal potential differences in cellular abundance in CTCL, we performed differential abundance analysis comparing CTCL with healthy skin, AD and psoriasis using milopy, a python implementation of Milo27. For the overall integrated object, we first performed a random subsampling using the subsample function in Scanpy, which subsampled the overall object to 0.2 of the total number of cells. Then the standard milopy pipeline was run for the data object with the number of neighbors (n_neighbors) being set to 100 and the proportion of graph vertices to randomly sample (prop in the milo.make_nhoods function) being set to 0.05. For the APCs, stromal and benign cell populations, we did not subsample the objects and set n_neighbors and prop to 100 and 0.05, respectively. Beeswarm plots were made using the plotDAbeeswarm function to show the log(transformed fold-changes) in abundance of cells in CTCL versus those in healthy skin, AD and psoriasis for each data object.

Inferring CNVs based on scRNA-seq data

To effectively distinguish malignant T cells and nonmalignant cells, we inferred large-scale chromosomal CNVs of single cells based on scRNA-seq data using the tool InferCNV (https://github.com/broadinstitute/inferCNV) with default parameters. In brief, InferCNV first orders genes according to their genomic positions (first from chromosome 1 to chromosome X and then by gene start position) and then uses a previously described sliding-average strategy to normalize gene expression levels in genomic windows with a fixed length. Multiple putative nonmalignant cells are chosen as the reference to further denoise the CNV result. From the CNV inference, we did not detect malignant T cells in six patients, which is consistent with the report from a study that we integrated21.

Intratumor expression programs and metaprograms

To explore intratumor expression programs, we applied non-negative factorization (implemented in the R NMF package) to the tumor cells from each patient. In brief, for each tumor, we first normalized the expression counts using the NormalizeData function in Seurat with default parameter settings. Highly variable genes (HVGs) were then selected using the FindVariableFeatures function in Seurat69. Next, we performed center scale for HVGs and regressed out the percentage of mitochondrial genes using the ScaleData function. For NMF analysis, all negative values in the expression matrix were replaced by 0. The top ten ranked NMF gene modules in each tumor sample were extracted using the nmf function in the NMF package. For each gene module, we extracted the top 30 genes with the highest weight, which were used to define a specific intratumor expression program. Finally, we included only intratumor expression programs that had s.d. values >0.1 among tumor cells. To investigate whether some intratumor expression programs were actually shared by multiple tumors, we applied a clustering analysis to all programs based on the pair-wised Jaccard index calculated as follows, where A and B represent two intratumor programs:

$$}=}\bigcap }/}\bigcup }.$$

We defined those intratumor programs shared by multiple tumors as metaprograms. We interpreted metaprograms by checking functions of individual genes and based on the 3CA database70 (https://www.weizmann.ac.il/sites/3CA).

Subclustering and annotation of different cell compartments

We performed subclustering and annotation of different cell compartments based on the object integrated with the skin cell atlas data. Cells were subsetted into stromal and keratinocyte, APC and benign T cell, NK cell and innate lymphoid cell populations. To improve separation of subclusters for the stromal and keratinocyte and APC populations, we reintegrated the cells to new embeddings. The top 2,000 and 1,000 HVGs were selected, respectively. To create new embeddings, an scVI model was set up with the donor as batch key, technology as technical covariate and the default scVI model parameters. We used the labels transferred from the skin cell atlas using CellTypist to annotate cells. These objects were then used to analyze differential abundance of granular cell types between CTCL and other conditions.

For the B cell compartment, as most B cells are from CTCL samples, we subsetted them from the CTCL object. HVGs were detected using the highly_variable_genes function in scanpy with minimum cut-off values of 0.0125 and 0.5 for expression and dispersion. We excluded immunoglobulin genes from HVGs when running principal component analysis. To correct the batch effect, we used the harmonypy71 package and set the donor as the batch key with the θ value being set to 3.

Visium processing and spatial mapping with cell2location

Sequencing reads from 10x Genomics Visium FFPE libraries were aligned to the human transcriptome reference GRCh38-2020-A using 10x Genomics SpaceRanger (v.2.1.0) and exonic reads were used to produce messenger RNA count matrices for each sample. The 10x Genomics SpaceRanger was also used to align paired histology images with mRNA capture spot positions in the Visium slide. To spatially map the cell types annotated in scRNA-seq data to their spatial locations in tissues, we applied cell2location to integrate scRNA-seq CTCL data and Visium FFPE mRNA count matrices as described previously28. Briefly, the cell2location model estimates the abundance of each cell type in each location by decomposing mRNA counts in Visium FFPE data using the transcriptional signatures of reference cell types derived from scRNA-seq data. Two major steps were in the analysis using cell2location: (1) we applied a negative binomial regression model implemented in cell2location and estimated the reference signature of fine-grained annotated cell types in integrated scRNA-seq data from CTCL and healthy samples. In this step, we used an unnormalized mRNA count matrix as input and filtered it to 12,711 genes and 626,529 cells. Donor IDs were regarded as the batch category, donor, study and 10x version as covariates, and the following parameters were used to train the model: ‘max_epochs’ = 250, ‘batch_size’ = 2,500, ‘train_ size’ = 1 and ‘Ir’ = 0.002. (2) The reference signature model was further used by cell2location to estimate spatial abundance of cell types. We kept genes that were shared with scRNA-seq and estimated the abundance of cell types in the Visium data from 8 CTCL samples (from 8 patients) and 15 healthy skin samples (from 7 donors). In this step, cell2location was used with the following parameter settings: training iterations = 50,000, number of cells per location N = 30, ‘detection_alpha’ = 20. We further performed NMF of the cell-type abundance estimates from cell2location to identify the spatial co-occurrence of cell types with ‘n_fact’ being set to 5–30. For downstream analysis, considering that our Visium samples consisted of a single piece of tissue, we removed all spots that correspond to tissue debris (which are all spot groups except the largest one). We considered outermost epidermis spots as skin surface and used them to calculate Euclidean distances from each spot to the closest surface spot and expressed this distance in interspot distances (100 μm). To estimate cell-type abundance dependent on distance to surface, we first normalized cell-type abundance by dividing it by per-spot totals. Then, we grouped spots by rounded distance to surface and calculated the mean and s.d. for each cell type and each distance.

DEG analysis using a pseudo-bulk strategy

We applied a pseudo-bulk strategy to the analysis of DEGs across (1) malignant and benign T cells, (2) malignant T cells from epidermis and dermis, (3) malignant T cells from early stage and advanced-stage samples, and (4) microenvironmental cells from CTCL and the other three conditions. In brief, we aggregated raw counts of each gene by donor and used donors rather than cells as biological replicates. DEG analyses were carried out using R package edgeR72. For the analysis of (1), we excluded nonlesion cells from AD and psoriasis and regarded healthy skin, AD and psoriasis as one comparator (other). We filtered genes by expression levels using the filterByExpr function in edgeR with ‘min.count’ and ‘min.total.count’ being set to default. We designed the model matrix using the model.matrix function and included only one variable, namely groups (malignant T cell and benign T cell). For the analysis of (2), we first divided malignant T cells from the eight patients with separate sampling of epidermis and dermis into those from epidermis and dermis, and conducted pseudo-bulk analysis with both tissue layer and derived patient being considered (for example, CTCL1_dermis and CTCL1_epidermis were aggregated separately). We filtered genes by expression levels using the filterByExpr function in edgeR with ‘min.count’ and ‘min.total.count’ being set to 2 and 10, respectively. In the model matrix design, we fit the model on paired samples considering both tissue (dermis and epidermis) and patient (CTCL1–CTCL8). For the analysis of (3), we included both dataset (Peking University (PKU), MD Anderson Cancer Center (MDA), Vienna, Ncl_Sanger_Fresh and Ncl_Sanger_FFPE) and group (early stage and late stage) as variables to consider variation across studies. For the analysis of (4), we included 10× CTCL data from fresh samples (both newly generated and published) because the probe-based 10x Flex for FFPE samples does not capture MHC-II genes. We excluded nonlesion cells from AD and psoriasis and regarded healthy skin, AD and psoriasis as one comparator (other). We excluded samples with <20 cells of a specific cell type. We filtered genes by expression levels using the filterByExpr function in edgeR, with ‘min.count’ and ‘min.total.count’ being set to default. We designed the model matrix using the model.matrix function and included only one variable, namely groups CTCL and other (healthy skin, AD and psoriasis). For all the analysis, we fit genewise, negative binomial, generalized linear models with quasi-likelihood tests using the glmQLFit and glmQLFTest functions in edgeR.

Bulk deconvolution of cell types

For bulk deconvolution analysis, we first downloaded published bulk RNA-seq datasets of healthy skin, AD, psoriasis, and CTCL from the Gene Expression Omnibus (GEO) database with accessions GSE121212 and GSE168508. A single-cell reference for deconvolution analysis was then prepared by excluding tumor cells and randomly downsampling the integrated object (healthy skin, AD, psoriasis and CTCL) to 8% of total cells. BayesPrism38 was used for deconvolution analysis with raw counts for both single-cell and bulk RNA-seq data as inputs. Both the ‘cell-type labels’ and the ‘cell-state labels’ were set to fine-grained annotations. Ribosomal protein genes and mitochondrial genes were removed from single-cell data because they are not informative in distinguishing cell types and can be a source of large spurious variance. We also excluded genes from sex chromosomes and lowly transcribed them as recommended by the BayesPrism tutorial. For further analysis, we applied a pairwise Student’s t-test to select DEGs with ‘pval.max’ being set to 0.01 and ‘lfc.min’ to 0.1. Finally, a prism object containing all data required for running BayesPrism was created using the new.prism() function, and the deconvolution was performed using the run.prism() function. Two-sided Wilcoxon’s rank-sum test was performed to examine any statistically significant enrichment.

For the survival analysis, the CTCL bulk RNA-seq cohort was grouped into high and low abundance of B cells (both estimated by bulk deconvolution and mean expression of CD79A and CD79B) by the optimal cut point determined using the cutp() function in the survMisc R package. We performed multivariate analyses using Cox’s proportional hazards model (coxph() function in the survival R package) to correct clinical covariates including age, gender and tumor stage for the survival analysis. Kaplan–Meier survival curves were plotted to show differences in survival time using the ggsurvplot() function in the survminer R package.

Inference of cell–cell interactions

We inferred potential cell–cell interactions using CellPhoneDB73 (v.4). In brief, we randomly downsampled the CTCL object to 100 cells per fine-grained cell type per donor. The generated object was then used to run CellPhoneDB analysis with default parameters and thresholds. For the downstream visualization, we used the plot_cpdb3 function in the R package ktplots (https://github.com/zktuong/ktplots).

Prediction of druggable targets using drug2cell

To predict potential druggable targets on B cells and malignant T cells, we ran drug2cell51 on these two cell types together with benign T cells as a comparator. Drug2cell is druggable target prediction tool that integrates drug–target interactions from the ChEMBL database (https://www.ebi.ac.uk/chembl) with single-cell data to comprehensively evaluate drug target expression in single cells. We first calculated per-cell scores of ChEMBL drug targets using the d2c.score() function. Then, we performed differentially expressed analysis on ChEMBL drugs by comparing B cells, benign T cells and malignant T cells using the scanpy tl.rank_genes_groups() function. When visualizing the result, we separated malignant T cells by patients to show drugs that potentially function in multiple patients, given the strong inter-patient heterogeneity of CTCL tumors.

IHC of FFPE samples for TOX and GTSF1

IHC staining for TOX and GTSF1 was performed on skin samples from healthy skin, AD, psoriasis and CTCL. In addition to the skin samples collected for scRNA-seq, a further cohort of patients with CTCL gave informed written consent for previous clinical samples to be used. Automated IHC staining was performed by NovoPath on FFPE slides using the Ventana Discovery Ultra autostainer (Roche) and the DISCOVERY ChromoMap DAB Kit. Antibodies used for staining were anti-GTSF1 (polyclonal, Atlas Antibodies, HPA038877) and anti-TOX (polyclonal, Atlas Antibodies, HPA018322). Scoring for TOX and GTSF1 was performed manually by a hematopathologist and dermatologist, reviewing the slides and deciding on an agreed approximation of positive staining. Identification of neoplastic T cells was based on their location, size and immunophenotype.

IHC of FFPE samples for B cells

IHC staining was performed on FFPE skin samples of CTCL tumors. Skin biopsies were fixed in 10% neutral buffered formalin (Newcastle and Stockholm cohorts) or 4% formalin (Vienna cohort), then moved to 70% ethanol, dehydrated and embedded in paraffin. For tissues from Vienna, FFPE samples were cut into 4-μm sections, deparaffinized using a Neoclear (Sigma-Aldrich) and ethanol series and autoclaved in citrate buffer at pH 6.1 (Dako) to achieve antigen retrieval. Blocking with hydrogen peroxide was performed. Subsequently, slide tissue sections were subjected to automated IHC staining (Autostainer, Dako Agilent) using anti-CD20 antibody (clone L26, Dako, M0755) or anti-CD79a antibody (clone JCB117, Dako, M7050) followed by visualization (EnVision FLEX, Dako Omnis, Agilent). For samples from Newcastle, 4-µm FFPE sections were stained for CD20 or CD79a on the Ventana Benchmark Ultra automated staining platform with CC1 antigen retrieval and Ultraview-DAB detection. For analyses in the Stockholm cohort, FFPE samples were cut into 3.5-µm sections, deparaffinized and subjected to citrate buffer at pH 9 (Dako) to achieve antigen retrieval. Staining was manually performed using anti-CD20 antibody (clone EP459Y, Abcam, ab78237) and secondary goat anti-rabbit immunoglobulin G H&L antibody (Abcam, ab214880).

Stained sections were imaged and digitized with NanoZoomer S360 (Hamamatsu), Scanscope CD2 (Aperio Technologies), Zeiss AxioScan.Z1 Slide Scanner and TissueFAXS scanning system (TissueGnostics). Image-based automated cell detection for all samples was performed with HistoQUEST software (TissueGnostics).

Statistics and reproducibility

For representative images from IHC and RareCyte, each experiment was repeated two to three times independently with similar results. Information on statistical tests is included in the figure legends.

Reporting summary

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