Characterisation of tumour immune microenvironment remodelling following oncogene inhibition in preclinical studies using an optimised imaging mass cytometry workflow

Mouse models are critical in pre-clinical studies of cancer therapy, allowing dissection of mechanisms through chemical and genetic manipulations that are not feasible in the clinical setting. In studies of the tumour microenvironment (TME), novel highly multiplexed imaging methods can provide a rich source of information. However, the application of such technologies in mouse tissues is still in its infancy. Here we present a workflow for studying the TME using imaging mass cytometry with a panel of 27 antibodies on frozen mouse tissues. We optimised and validated image segmentation strategies and automated the process in a Nextflow-based pipeline (imcyto) that is scalable and portable, allowing for parallelised segmentation of large multi-image datasets. Incorporating user-specific plugins, imcyto can be flexibly tailored to a wide range of segmentation needs. With these methods we interrogated the dramatic remodelling of the TME induced by a KRAS G12C inhibitor in an immune competent mouse orthotopic lung cancer model, showcasing their potential as key discovery tools to enhance understanding of the interplay between tumour, stroma, and immune cells in the spatial context of the tissue. presented here a complete IMC workflow for use in preclinical mouse studies and demonstrated the value and importance of using IMC to study the effects of a treatment on the spatial organisation of the TME.


Introduction
The advent of successful immune checkpoint inhibitors has revolutionised the treatment of cancer in recent years. However, a large proportion of patients exhibit intrinsic or acquired resistance, and good prognostic markers for response are lacking. The TME is thought to play a key role mediating immune evasion.
Therefore, it will be crucial to enhance our knowledge about the cells infiltrating the tumour and their spatial context. Manipulating the TME to revert immune suppression has the potential to significantly enhance the efficacy of immunotherapies. Mouse preclinical cancer models provide an excellent platform to study such interventions aimed at the TME in a controlled manner.
In recent years, with the use of multiplex techniques such as single cell mass cytometry (CyTOF) and single cell RNA sequencing, it has become apparent that tumours are infiltrated with a diverse spectrum of immune cells, often with different phenotypes from their normally homeostatic counterparts 1, 2 . Unfortunately, the digestion of the tissue that is required to perform such analysis destroys the spatial context of the TME. Immunofluorescence (IF) and immunohistochemistry can provide spatial localisation data, but the number of markers that can be used simultaneously is limited by the spectral overlap of fluorophores and chromogens.
Thus, novel highly multiplexed imaging technologies, such as imaging mass cytometry (IMC) that is based on unique atomic mass not wavelength, are extremely attractive, allowing for in depth characterisation of the TME with a metal-conjugated antibody panel of up to 40 markers while retaining the spatial context 3,4 .
While products and methods for conducting IMC studies on patient samples are becoming well established, publications using IMC in the mouse are scarce and use limited antibody complexity without subsequent quantification of the images 5, 6 . Here we present a complete IMC workflow, including a validated 27-marker antibody panel, automated and optimised image segmentation using our imcyto pipeline, and showcasing various spatial analyses. We applied these methods to study the effects of MRTX1257, a mutant specific inhibitor of the KRAS G12C oncoprotein, on the TME of an immunotherapy refractory KRAS G12C mutant lung cancer. The targeted inhibition of oncogenic KRAS signalling using the recently developed mutant specific van Maldegem, Valand et al. 5 March 2021 4 KRAS G12C covalent inhibitors has shown very promising efficacy in reports of currently ongoing phase 1 and phase 2 clinical trials in KRAS G12C mutant lung cancer 7,8 . However, it is expected that long term therapeutic responses will be limited by acquisition of drug resistance, and therefore combination therapies are under intense investigation. A recent study reported enhanced survival when combining the KRAS G12C inhibitor AMG510 with anti-PD1 immune checkpoint blockade in an immunotherapy sensitive syngeneic KRAS G12C mutant CT26 colon carcinoma subcutaneous tumour model 9 . This prompted us to investigate the effects of tumour-specific KRAS inhibition on the TME in the context of a preclinical model of lung cancer.
Using IMC, we have gained in-depth and quantitative insight into the phenotypes and spatial relationships of immune cells and stromal compartments of the mouse lung TME, and how KRAS G12C inhibition promotes remodelling into a more immune activated state.

An antibody panel for imaging mass cytometry on mouse tissues
Due to the limited availability of antibodies that are well validated for use on mouse formalin-fixed and paraffin-embedded (FFPE) tissues, especially with respect to immunology markers, we designed and validated a 27-antibody panel for use on frozen tissue sections. Example composite images are shown in Figure 1a and Supplemental Figure 1.
With this panel we are able to distinguish many different immune cell types that are thought to play a role in the TME, such as lymphocytes and various subsets of myeloid cells. In addition, this panel included markers to visualise the context of the tissue architecture, e.g. endothelium and fibroblasts, as well as phenotypic markers that describe maturation and activation state of both tumour and immune cells (Figure 1b and Supplemental Table 1). Twelve metals remain free for insertion of additional markers to customise this panel, with several slots in the most sensitive range of the detector, suitable for dimmer markers such as additional checkpoint molecules.

Targeted inhibition of KRAS G12C in a preclinical model of lung cancer
Orthotopic transplants of the 3LL DNRAS cell line, a KRAS G12C mutant and NRASknockout Lewis Lung Carcinoma derivative that we previously shown to be sensitive to KRAS G12C inhibition 10 , were treated for seven days with the KRAS G12C inhibitor MRTX1257 or vehicle, before lungs were harvested, processed and stained with our antibody panel. In marked contrast to CT26, the Lewis Lung Carcinoma is considered to be highly refractory to immune interventions 11 and has an immune cold, or immune excluded, phenotype. Treatment of 3LL DNRAS tumours established in the lungs with the KRAS G12C inhibitor MRTX1257 very markedly limited tumour growth over a period of seven days, compared to control mice, although most tumours did not actually reduce in size (Supplemental Figure 2).
A total of twelve tumours were selected for IMC, six tumours from three mice of each treatment group. At first inspection, the images revealed recurrent patterns in the cellular arrangement of the tissues (Figure 1c). Vehicle treated tumours usually van Maldegem, Valand et al.

March 2021
6 showed a tissue architecture in which CD68 + macrophages accumulated along the edge of the tumour, while F4/80 + and CD206 + macrophages seemed more intermixed with tumour cells. Effector cells such as CD4 + and CD8 + T cells were mostly excluded from the tumour area. Upon treatment with MRTX1257, the tissue organisation became more diffuse, macrophage phenotypes changed, and, most dramatically, the T cells infiltrated the tumour bed.

Optimisation of single cell segmentation using cell type specific settings
Quantification of IMC data requires image segmentation to extract expression data at the single cell level. As our dataset consisted of twelve large regions of interest (ROIs), ranging from 1 to 9 mm 2 , we sought a software solution that would be able to scale sufficiently and therefore chose to work with open source packages CellProfiler and Ilastik, as previously described by the Bodenmiller lab 12 . Any of the published segmentation strategies for IMC were originally optimised and validated on human FFPE tissues. Therefore, we first questioned what method would be most appropriate for use on frozen mouse tumours. In most instances, cell segmentation starts with the segmentation of the nuclei, which can subsequently be used as seeds to expand into cell objects. Various methods have been described to subsequently define the cell border: by expanding the nuclear perimeter by a fixed distance, such as 3 pixels, by a watershed method onto a membrane staining 4 , or by propagation onto a probability map of membranes generated by pixel classification with machine learning 12,13 . We started with testing two of these methods, the 3px expansion and the propagation method onto a probability map of membrane ( Figure 2a). From visual inspection of the cell masks we observed several recurring problems such as the failure to capture membrane staining of large cells, and the inclusion of signal from neighbouring cells into the smaller cells such as lymphocytes. In particular the inclusion of signal from neighbouring cells (spatial "spillover") led to increased noise in the data. In addition we observed that cells with a spindle shape, such as alpha smooth muscle actin (αSMA) expressing fibroblasts, were captured very poorly. We explored modifications of these segmentation strategies with the aim to improve the quality of the data. A simple strategy to improve the signal-to-noise ratio in the These three strategies used only one set of parameters to identify all nuclei, the seeds of each cell object. However, in complex tissues such as the TME, nuclei differ widely in size and shape. Setting a size threshold for an average tumour cell may easily exclude the smaller T cells, while stringent "declumping" settings can unnecessarily split up larger multi-lobed nuclei of tumour cells and macrophages.
Also, the general absence of a clearly identifiable nucleus associated with the αSMA signal was a likely cause for the poor segmentation of fibroblasts. Therefore, we modified the propagation method, to sequentially segment the different cell types using settings that are optimal for those cells. Each of the 'segmentation layers' is based on its own probability map and can use different settings for nucleus identification. Figure 2c is a graphic explaining the different layers of the sequential segmentation pipeline. By subtracting the mask of the cells identified at each step from the primary object layer (the nuclei), we can prioritise the cell types in the order at which they are segmented. The prioritised segmentation of the smaller lymphocytes improved their signal enrichment and signal-to-noise ratios, compared to the standard propagation method. A significant improvement was also seen for the capture of the αSMA signal into segmented cells (Figure 2b, d and e).
As a further validation of segmentation quality, we measured how well each method was able to detect challenging cell types, such as small lymphocytes, irregularly shaped macrophages, and the before mentioned fibroblasts. We manually annotated a set of CD4 + and CD8 + T cells, CD103 + dendritic cells (DCs), F4/80 + macrophages and αSMA + fibroblasts in the image. We counted how many of the highest expressing cells identified by the different segmentation strategies overlapped with these manually annotated cells. There was an improvement in detection of T cells, macrophages and fibroblasts in the sequential segmentation dataset, as shown in van Maldegem, Valand et al.

March 2021
8 the bar graph manual annotation in Figure 2f. As a negative control we also looked at the identification of CD103 + DCs, which was not prioritised in the sequential segmentation pipeline using cell type specific markers, but instead was segmented at step 9 along with the bulk of the remaining cells using 1px expansion. These cells were detected to a similar extent in all four segmentation pipelines.
An added advantage of using pixel classification with Ilastik was that we could also use marker combinations to identify larger tissue structures, including tumour, normal adjacent tissue and the interface between those two. This domain segmentation added a layer of spatial context, allowing for quantification of markers and cell types within the different tissue compartments (Figure 2g).

imcyto: A Nextflow based automated pipeline for segmentation on a computing cluster
To be able to run our complex sequential segmentation method on large datasets, we developed a Nextflow-based segmentation pipeline, optimised to run on a high performance computing cluster. Nextflow is based on the concept that software does not need to be installed on a server for it to be executed. Using Docker or Singularity containers, software packages are run as virtual images and are thereby portable, reproducible and platform independent. The resulting pipelines can be parallelised to make the process scalable across large datasets 14 . For automation of our workflow, we combined IMCtools, Ilastik and CellProfiler into our Nextflow-based "imcyto" pipeline (Supplemental Figure 3). The programs take in custom settings provided in the form of user-generated Ilastik and CellProfiler project files as plugins, thereby making the pipeline adaptable to widely varying segmentation methods. As such, all of the four segmentation strategies discussed above can be run using this pipeline.
Segmentation using the sequential method on our dataset consisting of 12 images ranging from 429MB up to 3.35GB, was completed in less than 45 minutes on the high performance cluster at The Francis Crick Institute.

IMC single cell analysis reveals dramatic remodelling of TME
The MRTX1257 treatment had clearly resulted in changes of phenotype and spatial organisation of cell types (Figure 1c). To interrogate the data in a more quantitative manner, we employed a collection of analysis methods, such as clustering and dimensionality reduction to characterise the cell types and their activation states. We Using CellProfiler's "Object relationships" output from our imcyto pipeline, we applied the Bodenmiller neighbouRhood analysis 15 on the two macrophage subtypes to look for cell types that would be enriched or depleted in their direct neighbourhood.
Rather than looking at the number of images for which a spatial relationship was van Maldegem, Valand et al. 5 March 2021 11 found significant, which is the default output of this analysis, we looked at the log2fold change in neighbourhood enrichment compared to the permutated data for each ROI, to reflect the magnitude of enrichment or depletion of cell types in each other's neighbourhoods (Figure 4e and Supplemental Figure 6). This revealed some interesting relationships, such as the consistent co-existence of Type 1 macrophages with DCs and lymphocytes, while they were being rarely found in close contact with tumour cells. Indeed, visualising these cell types showed that CD68 + macrophages are found often in close proximity to T cells and CD103 + DCs, particularly at the edge of the tumour (Supplemental Figure 5b). This is also true of the minority of intra-tumoural CD68 + macrophages, which maintain neighbourhoods with T cells and dendritic cells. Another interesting spatial interaction was observed for the Type 2 macrophages and fibroblasts, where in particular the fibroblasts in the tumour domain were found close to Type 2 macrophages (Supplemental Figure 5c). This is in line with reports that cancer associated fibroblasts can recruit monocytes and differentiate them to M2-like tumour associated macrophages 20 .
Lymphocytes are thought to be the most important effector in the anti-tumour  Figure 6a, b and c). To quantify in a more absolute way how the neighbourhood of the T cells changes with treatment, we decided to look at distances between cell types and frequency of those interactions. We calculated the distance of each cell to the nearest CD8 + T cell using Pythagoras's theorem, which gives the spectrum of proximities between cell types and works across the whole image. This approach showed that some cell types are consistently found in proximity to CD8 + T cells, such as other lymphocytes, dendritic cells and Type 1 macrophages. Other cell types were found further away at baseline, but this distance was reduced when the samples were treated with MRTX1257, such as for tumour cells and Type 2 macrophages -in agreement with the relocation of T cells into the tumour domain (Figure 4h). This result suggests that CD8 + T cells will come under stronger influence of tumour cells and Type 2 macrophages with treatment, which might relate to their upregulation of PD-1.
To use a more unbiased approach to dissecting the differences between vehicle and MRTX1257 treated samples, we performed multiple linear regression analysis on the mean intensities of all markers in the samples (Supplemental Table 2). As expected, active mTOR signalling as detected by the pS6 antibody showed the strongest reduction in response to MRTX1257 treatment, mainly from the inhibition of signal transduction pathways in the tumour cell clusters. Markers that were upregulated with MRTX1257 treatment included the macrophage maturation markers CD206, MHC class II, and F4/80, as discussed above, as well as the fibroblast marker αSMA. Interestingly, vimentin was the second highest coefficient in the table.
Vimentin expression is the canonical marker of epithelial-to-mesenchymal transition (EMT), which has been suggested as one of the possible mechanisms of resistance to KRAS-inhibition 21,22 . In fact, mean expression of vimentin increased for all cell types in the dataset, yet most striking was the difference seen for the blood vessels.

Discussion
Multiplex imaging techniques such as imaging mass cytometry are beginning to take centre stage in studies of the TME 24 . Application of this technology to pre-clinical mouse studies, however, has been lacking. Here we have described the design, optimisation and validation of a complete IMC workflow for mouse tissues. Our panel of 27 antibodies, provides a good basis for immuno-oncology studies, identifying many of the cell types of interest in the TME. An additional 12 isotopes are still available, to allow for customisation in any area of research interest, such as additional checkpoint molecules or stromal markers, or for combination with RNAscope, as was previously shown for human FFPE tissues 25 .
Considering the potential differences in cell densities and quality of the frozen mouse tumour tissues compared to human FFPE, we set out to find an optimised segmentation method, critical for obtaining good quality single cell data. We tested modifications to previously published image segmentation methods and established two strategies with a better performance with respect to signal-to-noise ratios and cell identification. The 1px expansion segmentation can represent a simple strategy to obtain quick segmentation data from a new antibody panel, without the need for training the classifiers. More flexible and sensitive segmentation can be obtained with a sequential segmentation strategy, and in particular this method can provide a solution to cells that are more challenging to segment. Any of those segmentation methods can simply be applied to large datasets, using our automated image segmentation pipeline, which has been made available to the wider community via the nf-core platform (nf-core/imcyto). This pipeline is scalable, yet flexible, as it can be customised to work with any antibody panel or segmentation strategy.
Working with mouse tissues, we have the advantage that we can image whole tumours, which is generally not possible with human samples. Seeing the tumour as a whole provides exceptional insight into the tumour architecture, as we have seen here for the very distinct cell communities within the different areas of the tissue, such as the accumulation of effector cells and Type 1 macrophages (CD68 + and CD11c + ) at the tumour edge whilst they were being excluded from the main body of the tumour. While MRTX1257 is a tumour-specific inhibitor, acting only on the G12C mutant form of KRAS protein found in the tumour cells and not the wild type KRAS protein found in the cells of normal tissue, it had profound indirect effects on the TME. The study of Canon et al. used flow cytometry and immunohistochemistry to look at changes in the TME upon KRAS G12C inhibition 9 . The main changes they observed were increases in macrophage, DCs and T cells. Our data confirmed those observations and provided additional spatial phenotypic data. We have used several different approaches to explore the spatial relationships between cell types, within cell neighbourhoods, and across the tissue domains. We showed how the tumour bed becomes infiltrated with maturing macrophages and PD-1 + lymphocytes, changing the spatial interactions between these cells. Interestingly some cells seemed to consistently exist together, such as the Type 1 macrophages that remain within the cell neighbourhood of T cells and DCs, irrespective of whether they are in the interface or in the tumour domain. This suggests either that their recruitment is regulated by the same mechanism, or that there is a functional importance to their interaction. While the DCs in that cellular niche expressed co-stimulatory molecules such as CD86 and could therefore potentially activate the T cells, the Type 1 macrophages expressed high levels of PD-L1 and would in turn be able to repress the T cell activation. Clearly, such local interactions could be fundamental to the outcome of immunotherapeutic treatments. It also raises the question of the state of the local chemokine and cytokine milieu. While a study by Schulz, et al. demonstrated that IMC can be combined with transcript detection using an RNAscope based in situ hybridisation 25 , the detection level of cytokine mRNA is often a limiting factor. A multi-omics approach combining IMC with single cell approaches such as CITE-Seq 26 , or spatial transcriptomics 27, 28 would therefore have the potential to significantly enhance our insight into the mechanisms regulating cellular communities.
A major advantage of a multiplex technology is that it allows for discovery and hypothesis generation. The finding of a higher expression of vimentin across all cell types was unanticipated but could be explained by a greater mobility of cells in the remodelling process of the TME 29 . In addition, it led to the observation that there was an increased presence of endothelium and possibly pericytes, pointing towards normalisation of the tumour vasculature upon MRTX1257 treatment in this normally haemorrhagic and necrotic tumour model 30 . Altogether, our data shows that van Maldegem, Valand et al.

March 2021
16 treatment with a tumour-specific KRAS inhibitor can dramatically remodel the TME in favour of an anti-tumour immune response, even in an immune cold tumour model such as the 3LL Lewis Lung Carcinoma. Current investigations are aimed at elucidating the mechanism by which this conversion is mediated.
We have presented here a complete IMC workflow for use in preclinical mouse studies and demonstrated the value and importance of using IMC to study the effects of a treatment on the spatial organisation of the TME.

In vivo drug study
This work was performed under a UK Home office approved project license and in accordance with institutional welfare guidelines. 10 6 3LL DNRAS Lewis Lung Carcinoma cells 10 were injected in the tail vein of 9-11 week old C57BL/6 mice and allowed to establish for 3 weeks. The lung tumour burden was assessed using a Quantum GX2 microCT Imaging System (PerkinElmer) and treatment groups were randomised with stratification by tumour number and size. MRTX1257 was prepared by sonication in 10% Captisol® (Ligand) and 50 mM citrate buffer (pH 5.0) and administered daily at 50mg/kg by oral gavage (5µg/g) for 7 days. At day 8 mice were scanned again and sacrificed with a terminal overdose (0.1ml/10g body weight intraperitoneal) of a mixture of Pentobarbital (2% w/v) and Mepivacaine hydrochloride (8 mg/ml), followed by cervical dislocation.

Tumour volume measurements
Mice were scanned before start of treatment and at the end of treatment using a Quantum GX2 microCT Imaging System (PerkinElmer), while mice were held under isoflurane anaesthesia. Tumour volume were measured using AnalyzeDirect software. Tumours with a minimum starting size of 0.4mm 3 at the time of the first scan were included in the analysis.

Tissue processing
Three mice from both the vehicle and the MRTX1257 treated group were processed to be used in IMC. Dissected lungs with tumours were stored in 20% ice cold sucrose up to one hour before embedding in Tissue-Tek O.C.T. Compound (Sakura) and being frozen gently using an isopentane liquid nitrogen bath. Blocks were stored at -80C until further processing.

Antibody staining
Five µm thin tissue sections were cut and collected on SuperFrost Plus™ Adhesion slides (Thermo Scientific) in a cryostat and stored at -80°C. When required, slides were thawed for 3 hours, fixed for 10 min with Image-iT™ Fixative Solution  Table 1.

Image acquisition
IF slides were imaged with a Zeiss Upright 710 microscope with a 20x objective lens using Zen blue imaging software (Zeiss). IMC images were acquired using Hyperion Imaging Mass Cytometer. Each ROI was selected such that it would contain a whole tumour including adjacent normal tissue where possible, or, if required, they were cropped to contain a single tumour using the cropping script (see Code availability).
Twelve images obtained from six mice were selected for this study, ranging 1~9 mm 2 (429Mb up to 3.35Gb).
For figures in this publication with IF or IMC images we used Fiji ImageJ v2.0.0 to make composite images of selected channels. For visualisation purposes the images were processed with an outlier removal step, filtered using a median or gaussian filter (0.5px radius) and scaled to enhance contrast. For Figure 1c the lymphocyte images were filtered using a band pass (3-40px) Fast Fourier Transformation to enhance detection of the cells.

Segmentation pipeline
The following section describes the proposed segmentation pipelines carried out in CellProfiler v3.1.9, including custom modules by Bodenmiller (https://github.com/BodenmillerGroup/ImcPluginsCP) and Ilastik v1.3.3b1 (see Code availability, and Supplemental Figure 3). In brief, data obtained as .txt files were converted into stacks of individual .tiff files per antibody marker using the IMCtools package v1.0.5 (https://github.com/BodenmillerGroup/imctools) 12 . Images in the full stack path were minimally filtered using outlier removal and median filtering in CellProfiler (both using the custom "Smooth Multichannel" module). For nuclear segmentation, images for 191Ir and 193Ir were summed, histogram equalised and segmented using a propagation-based thresholding strategy in CellProfiler.
For the pixel expansion strategy, segmented nuclei were expanded by a set number of pixels (1 or 3) to identify whole cell objects (see Figure 2a).
For the propagation strategy, selected markers (CD44, CD45, PECAM, CD11c, CD3, CD4, CD68, CD8, F480, aSMA, B220, CD206) were merged, alongside a nuclei image and converted into an RBG composite tiff in CellProfiler. This composite was fed into Ilastik pixel classification workflow and classified into nuclei, membrane and background to generate probability maps. In CellProfiler both segmented nuclei (as described) and membrane probability maps were used for whole cell segmentation using a propagation-based thresholding strategy (see Figure 2a).
Similarly, for the sequential segmentation strategy, multiple images representing individual cell types and domains were created by merging together selected markers in CellProfiler (PECAM for normal, CD44 for tumour, EPCAM and aSMA for structural, CD3, CD4, CD8 and B220 for lymphocytes, and CD68, F480, CD206 and CD11c for macrophages). These 5 new images alongside a nuclei image were combined into a 6-channel image in CellProfiler. This was fed into Ilastik autocontext workflow and classified into tissue domains (nuclei, tumour, normal, structural) and cell types (lymphocytes, macrophages, fibroblasts) to generate probability maps. In CellProfiler, individual cell types (in order -lymphocytes, macrophages, fibroblasts, normal cells, tumour cells, remaining cells) were sequentially segmented using both segmented nuclei (as described but with customised size parameters) and corresponding probability maps, with the exception of fibroblasts which were identified only using the probability map. All cell types were added together to generate a total cell mask. Domain segmentation was performed using thresholding

Validation of segmentation strategies
All cells in the ROI "BRAC3438.6f_ROI1_t1_Vehicle" (in short: "02_Vehicle"), were ranked for expression of each of the markers depicted, and the top 500 highest expressing cells were selected to calculate the mean intensity for that marker within the top 500, and the remaining cells, as a measure of signal enrichment. Signal-tonoise ratios were calculated by taking the mean intensity of relevant marker in the top 500 cells, relative to the expression of "polluting" markers not expected to be expressed in these same cells.
We generated a small manually annotated dataset from the same image set as above, by recording the x and y coordinates of CD4 + (n=92) and CD8 + T cells (n=70), CD103 + DCs (n=65), F4/80 + macrophages (n=108) and aSMA + fibroblasts (n=88). Cells from the top 500 expressing CD4, CD8, CD103, F4/80 or aSMA were interrogated to determine how many cells of each matched up to cells in the annotated dataset, following the criterium that the cell centre had to fall within 5px of each other.

imcyto pipeline
An automated Nextflow based pipeline that sequentially pre-processes and segments imaging data and extracts single cell expression information (see Supplemental Figure 2). This pipeline combines Docker or Singularity containers for IMCtools, CellProfiler and Ilastik. Inputs for the pipeline are image data in the form of .mcd, .txt, or ome.tiff files, a .csv file with binary information on the channels to be included, the user-generated CellProfiler pipeline files (.cppipe) -custom CellProfiler modules need to be separately provided, and (pre-trained) Ilastik project files (.ilp).
Output files from each step in the pipeline, such as raw or pre-processed TIFFs, probability maps, masks, Object relationships and measurements in the form of a .csv file will be created in an organised results folder structure. Alongside, a report is generated on the performance of the pipeline, with respect to memory and processor usage and running time for each individual processing step.

Phenotypic analysis
R implementations of tSNE (RtSNE) 31 and uMAP 32 were used for dimensionality reduction. Pseudotime analysis was done with Slingshot 33 , using 26 as start cluster.

Spatial analysis
A permutation approach developed by the Bodenmiller lab was used to determine whether detected neighbour interactions between metaclusters occurred more frequently in the images than observed by randomization (https://github.com/BodenmillerGroup/neighbouRhood) 15 . Briefly, an 'objectTable' logging object, and cluster information, and 'Object relationships' listing each cell and all its identified neighbours within a 15-pixel distance, were sent through 5,000 rounds of permutation. A P-value was then calculated to determine whether the difference between average baseline and permutated neighbour interactions were

Data availability
Image files, pipeline input files, single-cell data, 3D plot and associated files are available at https://hdl.handle.net/10779/crick.c.5270621.
The code generated during this study for analysis of single cell IMC data following image segmentation was written for R v3. 6    b. Representative false colour images of F4/80 (green), CD4 (yellow) and aSMA (red) markers merged with a nuclear stain (blue) and cropped (top row). Cell mask outlines (white) overlaid onto these markers were generated by either 3px expansion strategy (2nd row) or propagation strategy (3rd row), 1px expansion strategy (4 th row) or sequential segmentation (bottom panel). Image processing for visualisation: outliers were removed and a median filter of 0.5px radius was applied in Fiji.
c. Schematic of the sequential segmentation pipeline. Each step uses custom size threshold settings for nucleus detection, as well as cell type specific markers to generate the probability maps for membranes in Ilastik. A propagation step as described in (a) is subsequently used to find cell boundaries in steps 1, 2, 4, 5, 7. At every level, the identified objects are subtracted from the total remaining nuclei. Any