images—which implements DeepCell’s Mesmer 52 (v.0.4.1) for automated cell segmentation based on TissueNet, a dataset for training segmenta- tion models with more than one million manually labelled cells. In brief, Steinbock was run using the Steinbock Docker container in Linux as provided by the authors. Raw .mcd files were preprocessed and filtered for hot pixels with a threshold of 50. Next, cells were segmented using Mesmer based on CD11c, CD14, CD45, CD68, CD20, CD8a, CD169 and CSK-III as cytoplasm markers and SOX10, DNAI and DNAII as nuclear markers. Single-cell data were extracted from segmented images by calculating the mean object intensities per channel, as well as spatial object properties (area, centroid, major axis length, minor axis length and eccentricity). Processing and normalization of IMC data Data were processed using an R-based workflow for multiplex imaging data, as described previously 49 . In brief, data were transformed using an inverse hyperbolic sine function (asinh function from the rrscale v.1.0 package). The dataset then underwent median cell size calcu - lation, with very small cells (diameter less than 10 pixels) excluded from further analysis. Shot noise was filtered mostly by performing all downstream analyses using median signal intensities per cell, with cells consisting, on average, of 80–90 pixels per μ m 2 —minimizing the effect of rare random positive pixel artefacts. As a quantitative metric of image quality, we generated signal-to-noise ratio (SNR) values for all markers in our panel (Extended Data Fig. 1b). SNR values were calcu- lated using a two-component Gaussian mixture model for each marker for each cell, calculated as SNR = I s / I n , in which I s represents the signal intensity (mean intensity of pixels with a true signal) and I n represents the noise intensity (mean intensity of pixels containing noise). Using the Otsu thresholding method, we identified pixels of foreground (signal) and background (noise), in which signal intensity is divided by noise intensity for each cell as described previously 49 . Most hallmark markers for cell-type identification (for example, CD4, CD8, CD20, CD68 and CD11c) showed relatively strong signal intensities and SNR values (>2). Transcription factors and specific cell-surface receptors showed lower staining intensities, yet all had positive SNR values (> 1.5). ROI cover- age was calculated and any ROIs with a cell coverage lower than 50% were inspected to identify potential selection errors. To ensure data integrity, batch correction was then applied to the data per ablation run, as well as Ki-67 expression, using the fastMNN method 53 (v.1.29). Cell-type annotation FastMNN-corrected data were clustered using the Rphenograph pack- age (v.0.1.1) 54 , with k = 40. The resulting clusters were then annotated using marker proteins to distinguish general lineages (B cells, CD4 + T cells, CD8 + T cells, T reg cells, tumour, endothelial cells, myeloid cells, fibroblasts and stroma). These lineages were subclustered to remove cells that were wrongly annotated and to further stratify the lineages into more defined cell phenotypes. We refrained from using mark- ers with relatively low signal intensity and SNR values (Extended Data Fig. 1b) for phenotyping (for example, BCL6 and PD-L1). To increase annotation accuracy, we limited definitions of key cell types to naive and non-naive T cells (based on CD3, CD4, CD8A, FOXP3, CD45RA, Ki-67, CTLA4, PD-1 and HLA-DR); DCs (CD11c hi CD68 lo ) separated on the basis of high or low HLA-DR levels; and macrophages (CD11c lo CD68 hi ) separated on the basis of high levels of CD163 or CD169. The phenotyp- ing of cells was validated by careful inspection of labels projected on cell masks in combination with marker expression on images, often resulting in several rounds of visual inspection and subsequent param- eter optimization (for example, marker-based gating of individual cell phenotypes using cytomapper). Cell–cell interaction and neighbourhood analysis Interaction graphs were calculated using ‘buildSpatialGraph’ from the imcRtools package v.0.1.8 (ref. 49), with the Delaunay triangulation method, resulting in a definition of the nearest neighbours for each cell in the IMC dataset (maximum distance parameter = 20). Per ROI image, we used ‘testInteractions’ to compute the average cell-type–cell-type interaction counts and compared this count against an empirical null distribution generated by permuting all cell labels, while maintaining the tissue structure. This way, every pairwise interaction was tested for statistically significant interaction (+1) or avoidance (−1) in every ROI. Sums of interact and avoid scores were visualized in a heat map as per- centages of ROIs with a particular label. ‘testInteractions’ was also used to compute the average interaction count for each cell phenotype pair across the entire dataset, including label permutation to enable statisti- cal testing. Next, we inferred cellular NBs using ‘aggregateNeighbors’ based on the Delauney interaction graph, in which cells were counted by their phenotypic label and separated by disease stage as well as by location. Data were then k -means-clustered ( k = 10), from which the resulting NBs were visualized using a heat map for each patient group (stage II or III; R or DF) and location (for example, paracortex, intratumoural) separately. Cell distances were calculated using ‘minD- istToCells’ from the imcRtools package (v.0.1.8) 49 . Next, dyads were defined as surrounding cells within 15 μm of target cells as measured from the centroid of cells. We then calculated cell distances to cells labelled as DC–CD8 dyads and quantified macrophage abundance near these dyads (within 15 μm or 30 μm), separated by patient group. Image visualization Composite images of selected channels were generated using Fiji ImageJ (v.1.54t) 55 . Images were processed to remove outlier pixels, filtered using a median or Gaussian filter (0.5-pixel radius) and scaled to enhance contrast as described previously 56 , or visualized using the plotPixels function in the cytomapper 57 package. GeoMx digital spatial profiling Digital spatial profiling (DSP) using the GeoMx Whole Transcriptome Atlas (WTA) (NanoString Technologies) was performed following the manufacturer’s instructions. In brief, 5- μ m-thick slices were cut from TMAs and stained for SYTO 13 (NanoString Technologies, 1:20), CD11c–Alexa Fluor 647 (Abcam, EP1347Y, 1:40), CD8a–Alexa Fluor 594 (BioLegend, C8/144B, 1:40) and CD68–Alexa Fluor 532 (Novus Bio- logicals, C68/684, 1:40), and incubated overnight with GeoMx WTA probes. Within the TDLN paracortex, regions in which CD11c + , CD68 + and CD8 + cells co-localized were selected. From these regions, WTA probes from CD11c + , CD8 + and CD68 + cells were collected separately and sequenced using an Illumina NextSeq2000. Paired-end clusters were generated, producing reads that were 26 bases in length. We chose to remove segments with less than 15% of all genes detected. The resulting counts were normalized using Q3 normalization and DEGs ( P < 0.05) were calculated using a linear mixed model (LMM) as implemented in the GeomxTools v.3.10.0 R package (https://doi.org/10.18129/B9.bioc. GeomxTools). Pathway enrichment analyses on the DEGs were done using Metascape 58 . GSEA We ran GSEA as implemented by the clusterProfiler 59 (v.4.3.3) R pack- age against custom gene signatures for DCs 21,60–62 , macrophages 22 and CD8 + T cells 33,34 . The GSEA function was implemented with default settings, using a Benjamini–Hochberg correction for multiple testing of the resulting P values. Analyses of publicly available transcriptome data TCGA transcriptome data for skin cutaneous melanoma (SKCM) were obtained from the UCSC Xena portal (https://tcga.xenahubs.net). To calculate the levels of PLA2GD2 in various types of melanoma sample, we used log 2 ( x + 1)-transformed RNA-Seq by Expectation Maximization (RSEM)-normalized counts and curated phenotype data as directly available from the UCSC Xena portal. To investigate which cytokines