Tissue samples for scRNAseq
Healthy lymph nodes, unaffected by neoplastic disease, were acquired from donation after circulatory death (DCD) adult donors at the time of removal of organs for transplantation by the Cambridge Biorepository for Translational Medicine (CBTM) with ethical approval (reference 15/EE/0152, East of England—Cambridge South Research Ethics Committee) and consent from donor families. Lymph node samples were collected from inguinal, mesenteric, and thoracic regions at the end of the organ donation procedure, within 1-2 hours of cessation of circulation under cold ischaemic conditions. Lymph node specimens were maintained in ice-cold 0.9% saline for transport. Tissue was dissociated and processed as previously described26. Lymph node data from (290b, 298c, 302c) are shared with James et al. and were processed with additional flow sorting as described26. Patients who donated diseased tissue for scRNAseq were enrolled in the ‘Investigating how childhood tumors and congenital disease develop’ study (NHS National Research Ethics Service reference 16/EE/0394).
Tissue samples for Nanostring and Multiplexed Immunofluorescence
The study was approved by a UK NHS Health Research Authority (HRA) Research Tissue Bank (CEPA Biobank, Newcastle upon Tyne Hospitals NHS Foundation Trust, reference number: 17/NE0070). The Research Tissue Bank released link-anonymized, formalin-fixed and paraffin embedded (FFPE) tissue, from patients with Classic Hodgkin lymphoma. All tissue samples were surplus to diagnostic requirements and were released in accordance with the terms of the ethical approval (ref: 17/NE0070).
Single cell RNA sequencing analysis
Single-cell RNA sequencing data from Hodgkin lymphoma specimens was acquired by direct data transfer from the Stiedl lab (University of British Columbia) as a raw counts matrix11.
Mapping and quantification of scRNAseq generated for this study (Wellcome Sanger Institute (WSI) dataset) was performed using the kallisto-bustools toolkit (0.24.4)27, against the GRCh38 human reference. The unfiltered count matrix was profiled using the defaultDrops function in the DropletUtils package (1.10.3)28, in R (4.0.4).
Quality control of all cell-containing droplets was performed using the scanpy python package (version 1.7.1)29. Across the concatenated dataset (WSI dataset and Aoki et al. data), cells with fewer than 500 distinct transcripts, fewer than 600 total counts, or greater than 10% counts originating from mitochondrial genes were filtered (Extended Data Fig. 1b). Genes with fewer than 10 counts across the dataset were filtered. The doublet detection tool scrublet30 was run on a per-channel basis with default settings, before calculating the proportion of doublet calls per cluster in an initial high-resolution per-channel clustering solution (Extended Data Fig. 1c) . Highly variable genes were called using the scanpy toolkit setting n_top_genes = 3000, flavor=”seurat_v3”. We performed data integration with batch correction (using donor as the batch key) and dimensionality reduction concurrently by training a single-cell variational inference model using the scvi-tools python package (version 0.9.0a0)12. A 15-dimensional latent representation was used as input to nearest neighbor graph construction, Uniform Manifold and Approximation (UMAP), and Leiden clustering in scanpy.
Marker genes were calculated using the tf-idf metric and genes were ranked by the tf-idf metric, and p values calculated by a hypergeometric test, corrected for multiple testing (Benjamini-Hochberg method)31. These marker genes were used as a guide to manual annotation of cell types.
Cell types were annotated according to canonical marker expression (Extended Data Fig. 1d), supplemented with labelling using the celltypist python package (0.1.15)32. Concordance between annotated celltype labels and predicted reference annotations was assessed by calculating the Jaccard distance between cell label strings (Extended Data Fig. 1e-f).
Cell type predictions between the King et al. B cell atlas33 was calculated by training a multilayer perceptron classifier on the reference data and calculating predictions on unseen test data, implemented in the scikit learn package (0.24.2)34.
Relative abundance analysis was performed using the MiloR package (0.99.1)18. The annotated dataset was represented as a k=30 nearest neighbor graph calculated in scVI latent space. This graph representation was partitioned into 22,573 overlapping neighborhoods, and cell counts per 10X channel were calculated per neighborhood. Differential abundance between disease and health was calculated using the testNhoods function, setting the design as ~centre + disease, and norm.method=”TMM”.
Ligand receptor analysis was performed using CellPhoneDB (version 2.1.7)22. The count matrix was normalized to 1e4 total counts before input to the python implementation of CellPhoneDB which was run with --iterations=1000, --threshold=0.1, and --result-precision=3.
Microarray data analysis
The GEO accession GSE39133 was downloaded using the GEOquery R package (version 2.60.0)35. We performed differential expression between GC and HRSC samples using limma (version 3.48.1)36 and derived genesets by retaining genes with an absolute log fold change >3 and adjusted p value <0.01.
Transcription factor regulon and network analysis
Transcription factor regulons data were acquired using the dorothea R package (version 1.3.3)37. Regulons with confidence levels “A”, “B”, and “C” were used to calculate transcription factor activities using the run_viper command.
We derived a network graph of transcription factors and targets using data from the OmnipathR R package (version 3.12)38, plotting edges between transcription factors and targets, and between targets using the ggraph R package (version 2.0.5).
Nanostring spatial transcriptomics analysis
Ten of the FFPE cHL tissue biopsies used in the multiplexed IF analysis were selected for in situ transcriptomics (Extended Data Table 1), using the Nanostring GeoMx platform (Nanostring, Seattle, WA). Biopsies were selected on the basis of year of fixation (2014 onwards), high quality multiplex staining, and successful RNAscope quality control (QC) tests (Advanced Cell Diagnostics, Inc., Newark, CA). Five test slides were prepared, each arrayed with 2 unique, whole tissue biopsy sections; one slide also included reactive lymph node tissue. For each test slide, 4 consecutive 4-μm-thick sections were taken from each tissue block, prepared according to the manufacturer’s instructions, and numbered 1-4: Slide 1, Hematoxylin & Eosin counterstain; Slide 2, CD45 and PAX5 immunostains (fluorescent IHC); Slide 3, CD274 (PD-L1) and PDCD1 (PD-1) - RNAscope probes, DNA counterstain; Slide 4, Oligonucleotide RNA probe hybridization.
Circular regions of interest (ROIs; uniform 300-μm diameter) were labeled upon the CD274/PDCD1 RNAscope slides, with reference to H&E and IHC slides for histomorphological context. PD-L1 was selected as the CD274 probe fluorescent labelling was more robust than identifying HRSC cells by CD30 probe in initial experiments. In all selected tumors, the HRSC-dense areas of tissue were known to be PD-L1+ and HRSC were confirmed by DAPI signal and nuclear morphology. For each tumor, 3 ROIs were selected in ‘PD-L1high’ areas and 3 ROIs in ‘PD-L1low’ areas of tissue, unless otherwise stated. Assessments of PD-L1 intensity were made by visual inspection of normalized CD274 probe fluorescent signal within each tumor. Sclerotic bands, germinal centers, and areas of low cellularity were avoided in cHL tissues, by referring to DNA counterstaining and paired IHC. In the control RN tissue, 3 ROIs were labelled in germinal centers and 3 ROIs in interfollicular areas.
RNA oligonucleotides from the ‘Cancer Transcriptome Atlas’ (n = 1812 probes), coupled to photocleavable oligonucleotide tags (‘barcodes’), were hybridized overnight. The following day, hybridized barcodes were released from tissue by UV exposure from each ROI sequentially, before being counted and sequenced on the Illumina HiSeq 2500 platform.
Nanostring data analysis
Quality control steps were performed on exported data to exclude outlier probes and then data were normalized against the 75th percentile of signal from each tumor. Following data QC, 1372 (75.7%) probes were retained for analysis.
The nanostring expression matrix (ROI-by-gene) was analysed by principal component analysis (PCA) using the scanpy python package, before constructing a shared-nearest neighbor graph representation. Briefly, we identified k=10 nearest neighbors in PCA space using the FNN R package, before calculating the Jaccard distance between neighbor vectors after the approach of Levine et al39. The resulting graph was clustered using the Leiden algorithm40 with default settings.
Deconvolution was performed by first generating a matrix of one-hot-encodings of marker genes per cell type, including HRSC. We then used the rlm function in the MASS R package to perform robust linear regression. The resulting coefficients were treated as fractions to derive counts by multiplying by the number of nuclei in each ROI. These cell count estimates were then used to calculate differential abundance estimates for each cluster using DESeq241.
Differential expression analysis between clusters was performed using limma (version 3.48.1)36. The top 10 differentially expressed genes per cluster were selected for plotting (Extended Data Fig. 3c)
Tissue samples for multiplexed immunofluorescence
Formalin-fixed, paraffin embedded (FFPE) tumor biopsies (all excised lymph nodes) were selected from the pathology archives of Newcastle upon Tyne Hospitals NHS Foundation Trust, Newcastle upon Tyne, UK, by members of the Novopath (formerly CEPA) biobank, with appropriate ethical approval (reference number: 17/NE0070). Limited, link-anonymized clinical data for study patients were collected by the Novopath biobank. Hematoxylin and eosin-stained (H&E) tissue sections and immunohistochemical tissue sections were reviewed by a Haematologist, together with the original pathology reports. Fifty-four cases were selected for the study, including biopsies from patients of all ages and histomorphological subtypes (Extended Data Table 1). Patient biopsies were selected on the basis of high-quality, whole lymph node resection biopsy tissue. Patients with coexisting malignant diagnoses were excluded from the study, even if these were not apparent in the tissue.
Multiplexed immunofluorescence
Multiplexed immunofluorescence (IF) was performed by sequentially immunostaining 4-µm-thick tissue sections from selected FFPE tumor biopsies, with primary antibodies, secondary reagents, and unique fluorochromes. Each of the two multiplexes comprises 6 primary antibodies, followed by a nuclear counterstain, 4’,6-diamidino-2-phenylindole, as per published protocols10,42,43. All immunostaining was performed on the automated Ventana Benchmark platform (Roche Diagnostics, Rotkreuz, Switzerland). Slides were then air dried, mounted with Prolong Diamond Anti-fade mounting medium (#P36965; Life Technologies) and stored at 4˚C in a light proof box until image acquisition. The target antigens, antibody clones, and dilutions for markers included in this report are listed in Extended Data Table 2.
Image acquisition
Test regions for multiplex IF analysis were initially identified for each tumor biopsy by reviewing the H&E and immunofluorescence tissue sections at 4x magnification. Three non-overlapping regions of interest were then acquired at 20x magnification for each multiplex panel. Each of these 3 test regions comprised 9 contiguous fields of view (FOVs), arranged as a 3x3 rectangular grid (measuring 2008µm x 1508µm in total). Therefore, there were a total of 27 FOVs per tumor, for each multiplex.
Regions were selected to exclude tissue artefacts, best represent the overall tissue, and to include CD30+ HRSCs. Regions were imaged and deconvoluted using the Vectra multispectral imaging platform (Vectra 3 and Inform 2.4.8, Akoya Biosciences, Marlborough, MA), using specific spectral libraries.
Computational image analysis
Cell detection and phenotyping: cell nuclei were segmented from the DAPI channel using inForm (version 2.4.8, Akoya Biosciences, Marlborough, MA; last accessed October 2021) with a fixed size cytoplasm region of ≤3µm grown around each cell nuclei. Cells were phenotyped based on the combinations of functional markers expressed using the phenotyping tool provided by inForm. Additional signals were manually thresholded for each case (PD-L1, TIM-3, IDO1, LAMP3).
Visualization of spatial density distributions: to visualize the spatial distribution of each phenotype, we applied a kernel density estimate to the spatial coordinates of each cell phenotype and rendered using contour plots. The kernel bandwidth set for each phenotype using Scott’s method44.
HRSC Neighborhood Analysis: To examine the immediate microenvironment surrounding HRSC cells, we determined a HRSC-niche neighborhood based on identifying cells whose centroid lay within a 25µm distance from a HRSC centroid to represent the microenvironmental niche (Extended data Fig. 4a). We apply a minimum distance of 10 µm to minimise the effects of stain spillover between neighboring cells45. We calculate the cell proportions of each cell phenotype within the HRSC-niche and for the entire image as a basis for comparison. To establish whether there are statistically significant differences across the dataset, Q-Q plots were used to identify the distributions as non-normal, and Wilcoxon signed-rank test applied to calculate statistical significance (Extended data Fig. 4g).
Nearest neighbor network analysis: We examined whether each cell phenotype tended to colocalize with or avoid each other cell phenotype. For a query phenotype, A, and a distant phenotype, B, we matched each cell of phenotype A to its nearest cell of type B and calculated the distance between each pair of cells (the nearest-neighbor distance). To establish a baseline distance which is independent of the density of an individual phenotype, we also randomly permuted the cell labels of type phenotype B so that we could estimate expected distance if the spatial local of cells of phenotype B were randomly distributed within the sample. We report the ratio of the observed nearest-neighbor distance to the expected nearest-neighbor distance for each cell to identify the distribution (Extended data Fig 4h). The median of this distribution is used as a summary statistic to visualise the pairwise relationships between all phenotypes, with a median nearest-neighbor distance ratio less than 1 indicating that cells tend to be placed closer together in the sample compared to a random distribution, and values greater than 1 indicating they tend to be further away (Extended Data Fig. 4h).