A simple open-source method for highly multiplexed imaging of single cells in tissues and tumours

Intratumoural heterogeneity strongly influences the development and progression of cancer as well as responsiveness and resistance to therapy. To improve our ability to measure and analyze such heterogeneity we have developed an open source method for fluorescence imaging of up to 60 protein antigens at subcellular resolution using formalin-fixed, paraffin-embedded (FFPE) tissue samples mounted on glass slides, the most widely used specimens for the diagnosis of cancer and other diseases. As described here, tissue-based cyclic immunofluorescence (t-CyCIF) creates high-dimensional imaging data through successive acquisition of four-color images and requires no specialized instruments or reagents. We apply t-CyCIF to 14 cancer and healthy tissue types and quantify the extent of cell to cell variability in signal transduction cascades, tumor antigens and stromal markers. By imaging immune cell lineage markers we enumerate classes of tumour-infiltrating lymphocytes (TILs) and their spatial relationships to the tumor microenvironment (TME). The simplicity and adaptability of t-CyCIF makes it a powerful method for pre-clinical and clinical research and a natural complement to single-cell genomics.


INTRODUCTION
Advances in DNA and RNA profiling have dramatically improved our understanding of oncogenesis and propelled the development of targeted anti-cancer drugs 1 . Sequence data are particularly useful when an oncogenic driver is both a drug target and a biomarker of drug response, BRAF V600E in melanoma 2 or BCR-ABL 3 in chronic myelogenous leukemia, for example. However, in the case of drugs that act through cell non-autonomous mechanisms, such as immune checkpoint inhibitors (ICIs), tumour-drug interactions must be studied in the context of a multi-cellular environment that includes both cancer and non-malignant stromal and infiltrating immune cells. Multiple studies have established that these aspects of the tumor microenvironment strongly influence the initiation, progression and metastasis of cancer 4 and the magnitude of responsiveness or resistance to therapy 5 .
Single cell transcript profiling provides a means to dissect the tumour ecosystems and quantify cell types and states 6 . However, single-cell sequencing usually requires disaggregation of tissues, resulting in loss of spatial context 6,7 . Tissue imaging preserves this context but as currently performed, has relatively low dimensionality, particularly in the case of slides carrying formalin-fixed, paraffinembedded (FFPE) tissue slices. These are among the most commonly acquired clinical samples and are routinely used to diagnose disease following staining with Haemotoxylin and Eosin (H&E). The potential for immunohistochemistry (IHC) of such samples to aid in diagnosis and prioritization of therapy is well established 8 but IHC is primarily a single channel method: imaging multiple antigens typically involves sequential tissue slices or harsh stripping protocols (although limited multiplexing is possible using IHC and bright-field imaging 9 ). Antibody detection by formation of a brown diaminobenzidine (DAB) or similar precipitate is also less quantitative than fluorescence 10 . The limitations of IHC are particularly acute in immuno-oncology 11 in which it is necessary to quantify multiple immune cell types (regulatory and cytotoxic T cells for example) in parallel with expression of multiple tumor antigens and immune checkpoints, such as PD-1/PD-L1 and CTLA-4.
A variety of multiplexed approaches to analyzing tissues that do not involve conventional microscopy have been developed with the goal of simultaneously assaying cell identity, state, and morphology [12][13][14][15][16] . For example, FISSEQ 17 enables genome-scale RNA profiling of tissues at single-cell resolution and multiplexed ion beam imaging (MIBI) and imaging mass cytometry achieve a high degree of multiplexing and excellent signal to noise ratios using metals as labels and mass spectrometry as a detection modality 12,18 . Despite the potential of these new methods, they require specialized instrumentation and reagents, one reason that the great majority of translational and clinical studies still rely on H&E or single-channel IHC staining. There remains a substantial need for highly multiplexed methods that (i) minimize the requirement for specialized instruments and costly, proprietary reagents, (ii) work with conventionally prepared FFPE tissue samples (iii) enable imaging of 20-60 antigens at subcellular resolution across a wide range of cell and tumour types (iv) collect data with sufficient throughput that large specimens (several square centimeters) can be imaged and analyzed (v) generate high resolution data typical of optical microscopy and (vi) allow investigators to customize the antibody mix to specific questions or tissue types. Among these requirements the last is particularly critical: at the current, early stage of development of high dimensional histology, it is essential that pathologists be able to test the widest possible range of antibodies and antigens in search of those with the greatest diagnostic and predictive value.
This paper describes an open-source method for highly multiplexed fluorescence imaging of tissues, tissue-based cyclic immunofluorescence (t-CyCIF), that significantly extends a method we previously described for tissue culture cells 19 . t-CyCIF assembles images of FFPE tissue slices stained with up to 60 different fluorescent antibodies via successive rounds of 4-channel imaging. t-CyCIF uses widely available reagents, conventional slide scanners, automated slide processors and freely available protocols to create a method that is easy to implement in any research or clinical laboratory. We believe that high dimensional imaging methods by t-CyCIF will become a powerful complement to single cell genomics, enabling routine analysis of the phenotypic geography of cancer at single-cell resolution.

t-CyCIF enables multiplexed imaging of FFPE tissue and tumor specimens at sub-cellular resolution.
In t-CyCIF, multiplexing is achieved using an iterative process (a cycle) performed on a slice cut from a block of FFPE tissue. We have developed t-CyCIF staining conditions for antibodies targeting ~140 different proteins including immune lineage makers, signaling proteins and phosphorylated kinases and transcription factors (some of which are drug targets), markers of cell state including cell cycle stage, quiescence and apoptosis etc. (Table 1). In the implementation described here, each cycle involves four sequential steps ( Figure 1A): (i) staining with fluorophore-conjugated antibodies against different protein antigens; we currently use antibodies conjugated to Alexa 488, 555 and 647 (ii) staining with Hoechst 33342 to mark nuclei (iii) four-channel imaging at low and high magnification (10X and 40X objectives) (iv) fluorophore oxidation using hydrogen peroxide, high pH and UV light followed by a wash step. To reduce the level of auto-fluorescence and minimize non-specific staining, we perform a pre-staining cycle prior to any incubation with primary antibodies ( Figure S1A); pre-staining involves incubation with secondary antibodies alone followed by fluorophore oxidation ( Figure 1B, Figure S1B-D). The current protocol has been optimized for samples prepared in the standard manner for pathologic diagnosis of cancer (4-5 µm thick FFPE slices mounted on a glass slide). We find that incubation in oxidation solution for 15 min is often adequate for Alexa 555 and 647 conjugated antibodies but reduction of Alexa 488 fluorescence to background levels typically requires 60 min ( Figure 1C-E, Figure S1E-G); we routinely perform 60-minute oxidation reactions.
Virtually all tissue samples we have examined can be successfully subjected to 8-20 t-CyCIF cycles, yielding data on the spatial distributions of 24-60 different antigens plus nuclear morphology.
The primary requirement appears to be good cellularity: samples in which cells are very sparse tend to be too fragile for repeated imaging, We achieve subcellular resolution using a fluorescence slide scanner (in our studies a RareCyte CyteFinder with 10X 0.3 NA objective and field of view of 1.6 x1.4 mm and a 40X 0.6 NA objective and field of view of 0.42 x 0.35 mm). Immunogenicity does not fall appreciably with cycle number ( Figure 1F-G) and the signal-to-noise ratio can actually increase due to a reduction in auto-fluorescence as cycle number increases ( Figure S1H-J) 20 . The primary limitation in the number of cycles is tissue integrity: some tissues are physically more robust and can undergo more staining and washing procedures than others ( Figure 1H). To date the highest cycle number has been obtained with normal tonsil and skin and with glioblastoma multiforme (GBM), pancreatic cancer and melanoma ( Figure S2 and S3). Cell morphology is preserved through multiple cycles: for example, in a t-CyCIF image of tonsil tissue ( Figure 2A, Table S1), we can distinguish membrane staining of anti-CD3 and CD8 in Cycle 2, and staining of the nuclear lamina, and nuclear exclusion of NF-ĸB in Cycle 6 ( Figure 2B).
Multiplexed immunofluorescence enables high resolution imaging of large samples. Figure 3 shows a ~2 x 1.5 cm pancreatic ductal adenocarcinoma (PDAC) subjected to 8 rounds of t-CyCIF ( Figure 3A, Table S2). The image comprises ~140 10X fields stitched together to reconstruct the full specimen. Differences in subcellular distribution are evident for many markers, but for simplicity, in this paper we only analyze intensity values integrated over a whole cell. Images were segmented using a conventional watershed algorithm and total fluorescent signal was calculated for each cell and antibody stain (Figure 3B-C; see Online Methods for a discussion of caveats). This yielded ~1.5 x 10 5 single cells each with 25 intensity values. When we analyzed the levels of pERK T202/Y204 (henceforth pERK, the phosphorylated, active form of the kinase) on a cell by cell basis we found that they were highly correlated with the levels of an activating phosphorylation of the downstream kinase pS6 S235/S236 (r = 0.81). Similarly, β-catenin levels (a measure of canonical WNT pathway signaling) were highly correlated with E-cadherin and Keratin levels, whereas Vimentin and VEGFR2 receptor levels were anti-correlated ( Figure 3D), recapitulating the known dichotomy between epithelial and mesenchymal cell states in normal and diseased tissues The WNT pathway is frequently activated in PDAC and is important for tumourigenesis of multiple gastrointestinal tumours 21 . Approximately 90% of sporadic PDACs also harbor driver mutations in KRAS, activating the MAPK pathway and promoting tumourigenesis 22 . Studies comparing these pathways have come to different conclusions with respect to their relationship: some studies show concordant activation of MAPK and WNT signaling and others argue for exclusive activation of one pathway or the other 23 . Using t-Stochastic Neighbor Embedding (t-SNE), which clusters cells in 2D based on their proximity in the 25-dimensional space of image intensity data, we identified multiple subpopulations within the same tumor sample representing negative, positive or no correlation between pERK and β-catenin levels (marked with labels "a", 'b" or "c", respectively in Figure 3E). In PDACs malignant cells can be distinguished from stromal cells, to a first approximation, by high proliferative index, which we assessed by we measured by staining for Ki-67 and PCNA 24 . When we gated for cells that were both Ki67 high and PCNA high cells and thus likely to be malignant, we again failed to find a fixed relationship between pERK and β-catenin levels. While we cannot exclude the possibility of phospho-epitope loss during sample preparation, it appears that the full range of possible relationships between the MAPK and WNT signaling pathways described in the literature can be found within a tumor from a single patient. This illustrates the impact of tumor heterogeneity on the activities of key signal transduction pathways.

Multiplex imaging of tumour architecture and immune infiltrates
Immuno-oncology drugs, including ICIs that target CTLA-4 and PD-1/PD-L1 are rapidly changing the therapeutic possibilities for traditionally difficult-to-treat cancers, such melanoma, renal and lung cancers, but responses are still highly variable across and within cancer types. Expression of PD-L1 correlates with responsiveness to ICIs such as pembrolizumab and nivolumab 25 but the negative predictive value of PD-L1 alone is insufficient to stratify patient populations 26 . Imaging for multiple markers such as PD-1, PD-L1, CD4 and CD8 using IHC on sequentially cut tumour slices appears to represent a superior approach to identifying ICI-responsive metastatic melanomas 5 . As a first step in developing multiplexed CyCIF immune biomarkers we developed staining conditions for antibodies against CD45, CD3, CD8, CD45RO, FOXP3, PD-1 and PD-L1 (Table 1). In a typical FFPE section from clear-cell renal cell carcinoma ( Figure 3F) we could distinguish a domain rich in non-malignant stroma, which stained strongly for the alpha isoform of smooth muscle actin (α-SMA) and one enriched in tumour cells and low in α-SMA. To examine the spatial distribution of cytotoxic CD8+ T cells within this tumor we stained for several immune markers, including CD8 ( Figure 3G), CD3, PD-1 and PD-L1. CD3 + CD8 + TILs were 4-fold enriched in the tumour-domain ( Figure 3H). The subset of CD3 + CD8 + PD-1 + cells that represents a population of putatively exhausted T cells was 18-fold enriched in the tumor-domain ( Figure 3H). Moreover, PD-1 and PD-L1 positive cells were 13 to 20-fold more prevalent in the tumour-rich domain (yellow bars, α-SMA low domain) as compared to the tumor stroma (blue bars, α-SMA high domain). Together, these data suggest that PD-1/PD-L1 interactions occur predominantly within tumour rich domains of kidney cancer and show the potential of t-CyCIF to quantify key features of TILs in combination with their positions in a tumour. Showing that such features constitute true biomarkers will, of course, required additional validation with clinical cohorts.

Analysis of diverse tumour types and grades using CyCIF of tissue-microarrays (TMA)
Tissue microarrays (TMAs) provide a means to analyze a large number of tissues and tumour samples simultaneously, and are frequently used to study patient biopsies collected in clinical trials. We applied 8 cycle t-CyCIF to TMAs containing 39 individual biopsies from 13 healthy tissues as well as low and high-grade tumours for 13 type of cancer ( Figure 4A, Figure S4, Table S2 for antibodies used, Table S3 for TMA details and naming conventions) and then performed t-SNE on single cell intensity data ( Figure 4B). The great majority of TMA samples mapped to one or a few discrete locations in the t-SNE projection (compare normal kidney tissue -KI1, low grade tumours -KI2, and high grade tumours -KI3; Figure 4C), while other tumours, such as ovarian cancer, showed a scattered pattern in the t-SNE projection ( Figure 4D). Overall, there was no separation between normal tissue and tumours regardless of grade ( Figure 4E). In a number of cases, high grade cancers from multiple different tissues of origin co-clustered, implying that transformed morphologies and cell states were closely related. For example, while healthy and low grade pancreatic and stomach cancer occupied distinct t-SNE domains, high grade pancreatic and stomach cancers were intermingled and could not be readily distinguished ( Figure   4F), recapitulating the known difficulty in distinguishing high grade gastrointestinal tumours of diverse origin by histophathology. 27 Nonetheless, t-CyCIF might represent a means to identify discriminating biomarkers by efficiently sorting through large numbers of alternative antigens, particularly those informed by known genetic features of each disease. Overall we conclude that t-CyCIF can be used on a wide range of normal tissues and tumor types to quantify inter-and intra-tumour heterogeneity.
Quantitative analysis reveals global and regional heterogeneity and multiple histologic subtypes within the same tumour in glioblastoma multiforme (GBM).
Data from single-cell genomics has revealed the extent of intra-tumour heterogeneity 28 but our understanding of this phenomenon would benefit greatly from spatially-resolved data 12 . To study this using t-CyCIF we performed 8 cycle imaging on glioblastoma multiforme, a highly aggressive and genetically heterogeneous 29 brain cancer that is classified into four histologic subtypes 30 . We imaged a 2.5 cm x 1.8 mm resected tumor sample for markers of neural development, cell cycle state and signal transduction state ( Figure 5A-B, Table S5). Phenotypic heterogeneity at the level of single tumor cells was assessed at three spatial scales corresponding to: (i) 1.6 x 1.4 mm fields of view (252 total) each of which comprised 10 3 to 10 4 cells (ii) seven macroscopic regions of ~10 4 to 10 5 cells each, corresponding roughly to tumour lobes and (iii) the whole tumour comprising ~10 6 cells. To quantify local heterogeneity we computed the informational entropy on a-per-channel basis for 10 3 randomly selected cells in each field ( Figure 5C, see online Methods for details). In this setting, informational entropy is a measure of cell-to-cell heterogeneity on a mesoscale corresponding to 10-30 cell diameters. For a marker such as EGFR, which is a driving oncogene in GBM, informational entropy was high in some areas ( Figure 5C; red dots) and low in others (blue dots). Areas with high entropy in EGFR abundance did not co-correlate with areas that were most variable with respect to a downstream signaling protein such as pERK. Thus, the extent of local heterogeneity varied with the region of the tumor and the marker being assayed.
Unsupervised clustering using expectation-maximization Gaussian mixture (EMGM) modeling on all cells in the tumour yielded eight distinct clusters, four of which in aggregate encompassed 85% of the cells ( Figure 6A). Among these, cluster one had high EGFR levels, cluster two had high NGFR and Ki67 levels and cluster six had high levels of vimentin; cluster five was characterized by high keratin and pERK levels. The presence of four highly populated t-CyCIF clusters is consistent with data from single-cell RNA-sequencing of ~400 cells from five GBMs 7 . Three of the t-CyCIF clusters have properties reminiscent of classical (cluster 1), pro-neural (cluster 2) and mesenchymal (cluster 6) histological subtypes, but additional work will be required to confirm such assignments. To study the relationship between phenotypic diversity and tumor architecture, we mapped each cell to an EMGM cluster (denoted by color). Extensive intermixing was observed at the level of fields of view and overall tumor domains ( Figure 6B). For example, field of view 147 was highly enriched for cells corresponding to cluster 5 (yellow), but a higher-magnification view revealed extensive intermixing of four other cluster types on a scale of ~3-5 cell diameters ( Figure 6C). At the level of larger, macroscopic regions, the fraction of cells from each cluster also varied dramatically ( Figure 6D, Figure S5). These findings have several implications. First, they suggest that GBM is a phenotypically heterogeneous on a spatial scale of 5-1000 cell diameters and that cells corresponding to distinct t-CyCIF clusters are often found in the vicinity of each other. Second, sampling a small region of a large tumour has the potential to substantially misrepresent the proportion and distribution of tumour subtypes, with implications for prognosis and therapy. Similar concepts likely apply to other tumor types with high genetic heterogeneity, such as metastatic melanoma, as recently indicated by single-cell genomic analyses 6 , and are therefore relevant to diagnostic and therapeutic challenges arising from tumor heterogeneity.

DISCUSSION.
The spatial heterogeneity of solid tumours poses a scientific, diagnostic and therapeutic challenge that is not sufficiently addressed using current methods. We have developed a simple, publicdomain approach for quantitative assessment of 20-60 protein antigens in ~5-10µm thick FFPE tissue slices, which represent the norm for diagnosis of human disease and study of mouse models. We describe several applications of t-CyCIF in studying oncogenic signaling, tumour heterogeneity and immune cell-tumour interaction, none of which requires specialized equipment (beyond a slide scanner) or proprietary reagents. t-CyCIF is not as technically sophisticated as FISSEQ 17 , MIBI 18 or tissue-based mass cytometry 12 , but we regard simplicity as a primary virtue: we have taught t-CyCIF to several other research groups and are confident that it can be readily adopted by many clinical and translational research laboratories. Cyclic immunofluorescence also appears to be substantially higher in throughput than non-optical methods, particularly when multiple slides are processed in parallel, and considerable opportunity exists for further improvement, for example, by switching from four to six channel imaging per cycle. Good linearity and resolution (~ 400 nm laterally in the current work, but potentially better with higher NA optics or super-resolution microscopes) are additional advantages of direct fluorescence imaging as compared to methods that rely on enzymatic amplification, laser ablation or mechanical picking of tissues. However, some signals, particularly those associated with phospho-epitopes can be very dim and amplification or use of very bright fluorophores such as quantum dots will be required to image them. We therefore continue to optimize t-CyCIF and will make updates available via our website.
In an initial test of t-CyCIF we quantified the relationship between WNT and MAPK-signaling Thus, what appears to be a set of conflicting findings most likely represents heterogeneity arising from differences in microenvironment, genotype or both. Such data adds new insight into our understanding of disease mechanism but variability may complicated the use of MEK-inhibitors in PDAC. 32 In a second test of t-CyCIF we studied within-tumor heterogeneity in GBM, a cancer with multiple histological subtypes whose differing properties impact prognosis and therapy. 30,33 Clustering antigen abundance data from t-CyCIF images also reveals multiple phenotypic classes within a single tumor. We have not yet established the link between t-CyCIF clusters and known histological subtypes but our results show that cells with very different characteristics are intermingled at multiple spatial scales with no evidence of recurrent patterns of heterogeneity. In regions of the GBM we have studied in detail, heterogeneity on a scale of 10-100 cell diameters is as great as it is between distinct lobes. The proportion of cells from different clusters also varies dramatically from one lobe to the next. Variation on this spatial scale is likely to impact the interpretation of small core needle biopsies.
In a third test we characterized tumour-immune cell interactions. Immune checkpoint inhibitors (ICI) produce durable responses in a portion of patients but identifying potential responders and nonresponders remains a challenge. Conventional single channel IHC on checkpoint proteins and ligands lacks sufficient positive and negative predictive value to stratify therapy or justify withholding ICIs in favor of small molecule therapy. Multivariate predictors based on multiple markers such as CD3, CD4, CD8, PD-1 etc. are likely to be more effective, but still underperform as patient stratification approaches 5 . This likely arises because cell types other than CD8+ TILs, including malignant, stromal, and myeloid-derived cells affect responses to ICIs. t-CyCIF represents a simple method for simultaneously assessing up to 60 predictors, including several processes and mechanisms, such as angiogenesis regulators, DNA damage, tyrosine kinase expression, proliferation, and others ( Figure S3).
In the current work, we perform a simple analysis showing that TILs can be subtyped, analyzed for PD-1 levels and proximity to PD-L1 ligand at a single cell level. Next steps include validating more antibodies and developing an efficient means to relate staining patterns to immune cell classes that have been defined primarily by flow cytometry.
Highly multiplexed histology is still in an early stage of development and better methods for segmenting cells, quantifying fluorescence intensities and analyzing the resulting data are required. With better data analysis methods, cell-to-cell heterogeneity in t-CyCIF images should enable reconstruction of signaling network topologies relevant to different regions of a single tumor 12,34 based on the observation that proteins naturally and randomly fluctuate in abundance and activity from one cell to the next. When these fluctuations are highly correlated, they are likely to reflect causal associations 35 .
Additional work will also be required to determine which types of heterogeneity are most significant for therapeutic response and disease progression: some cell-to-cell differences observed by fixed cell imaging arise from time-dependent fluctuation. However, we have observe cell-to-cell heterogeneity not only among proteins known to change over the course of a single cell cycle but also among long-lived proteins. Finally, validation of t-CyCIF-based biomarkers will require extensive testing in patient cohorts. It is important to note in this context that the current study describes a technology and its potential applications to histopathology, not actual diagnostic biomarkers. The associations we describe might not prove statistically significant when tested on larger, well-controlled sets of samples.
In conclusion, the t-CyCIF approach to multi-parametric imaging is robust, simple and applicable to many types of tumours and tissues; as an open platform, it allows investigators to mix and match antibodies depending on the requirements of a specific type of sample. In the long run we expect t-CyCIF to be complementary to, and used in parallel with more sophisticated protein and RNA imaging methods that may have greater sensitivity or channel capacity, although it seems probable that direct imaging will always have better resolution and speed than laser ablation or mechanical picking. A particularly important future development will be cross-referencing tumor cell types identified by singlecell genomics with those identified by multiplexed imaging, making it possible to precisely define the genetic geography of human cancer.

Competing financial interests
PKS is a member of the Board of Directors of RareCyte Inc., which manufactures the slide scanner used in this study, and co-founder of Glencoe Software, which contributes to and supports open-source OME/OMERO image informatics software. Other authors have no competing financial interests to disclose.   Table S1 for a list of antibodies used. (B) Images of selected channels and cycles emphasizing sub-cellular features.  Numbers refer to sample type; "1" to normal tissue, "2" to -grade tumors and "3" to high grade tumors.

ONLINE METHODS
We continue to improve the methods described here; periodic updates can be found at our web site http://lincs.hms.harvard.edu/resources/. A video illustrating the t-CyCIF approach can be found at https://youtu.be/fInnargF2fs.

Patients and specimens
Tumor tissue and FFPE specimens were collected from patients under IRB-approved protocols (DFCI

Reagents and antibodies
All conjugated and unconjugated primary antibodies used in this study are listed in Table 1

Automated dewaxing, rehydration and pre-staining
Pre-processing of FFPE tissue and tumor slices mounted on slides was performed on a Leica BOND RX automated stained using the following protocol: Step Steps 20-28 Pre-staining procedures as shown in Figure S1A: Step 20: IF Block -Immunofluorescence blocking in Odyssey blocking buffer (LI-COR, Cat.

927401)
Step 21: Antibody Mix -Incubation with secondary antibodies diluted in Odyssey blocking buffer Step 25: Staining with Hoechst 33342 at 2 μg/ml (w/v) in in Odyssey blocking buffer Manual dewaxing, rehydration and pre-staining. In our experience dewaxing, rehydration and prestaining can also be performed manually with similar results. For manual pre-processing, FFPE slides were first incubated in a 60°C oven for 30 minutes. To completely remove paraffin, slides were placed in a glass slide rack were then immediately immersed in Xylene in a glass staining dish (Wheaton 900200) for 5 min and subsequently transferred to a another dish containing fresh Xylene for 5 min.
Rehydration was achieved by sequentially immersing slides, for 3 min each, in staining dishes containing 100% ethanol, 90% ethanol, 70% ethanol, 50% ethanol, 30% ethanol, and then in two successive 1xPBS solutions. Following rehydration, slides were placed in a 1000 ml beaker filled with 500 ml citric acid, pH 6.0, for antigen retrieval. The beaker containing slides and citric acid buffer was microwaved at low power until the solution was at a boiling point and maintained at that temperature for 10 min. After cooling to room temperature, slides were washed 3 times with 1xPBS in vertical staining jars followed by blocking with Odyssey blocking buffer. Buffer was applied to slides as a 250-500 μl droplet for 30 mins at room temperature; evaporation was minimized by using moist in a slide moisture chamber (Scientific Device Laboratory, 197-BL). Slides were then pre-stained by incubation with diluted secondary antibodies (listed above) for 60 minutes, followed by washing 3 times with 1x PBS.
Finally, slides were incubated with Hoechst 33342 (2 μg/ml) in 250-500 μl Odyssey blocking buffer for 30 min. in a moisture chamber and washed 3 times with 1xPBS in vertical staining jars.

Cyclic immunofluorescence with primary antibodies and Hoechst 33342
All primary antibodies (fluorophore-conjugated and unconjugated) were diluted in Odyssey blocking buffer. Slides carrying dewaxed, pre-stained tissues were then stained in a moisture chamber by dropping the diluted primary or fluorophore-conjugated antibody (250-500 μl per slides) on tissue followed by incubation at 4 o C for ~12 hr. Slides were washed four times in 1x PBS by dipping in a series of vertical staining jars. For indirect immunofluorescence, slides were incubated in diluted secondary antibodies in a moisture chamber for 1 hr at room temperature followed by four washes with 1x PBS. Slides were incubated in Hoechst 33342 at 2 μg/ml in Odyssey blocking buffer for 15 min at room temperature, followed by four washes in 1x PBS. Stained slides were mounted prior to image acquisition (see the Mounting section below).

Primary antibodies
For t-CyCIF, we selected commercial antibodies previously validated by their manufacturers for use in immunofluorescence, immunocytochemistry or immunohistochemistry (IF, ICC or IHC). When possible, we checked antibodies on reference tissue known to express the target antigen, such as immune cells in tonsil tissue or tumor-specific markers in tissue microarrays. The staining patterns for antibodies with favorable signal-to-noise ratios were compared to those previously reported for that antigen by conventional antibodies. An updated list of all antibodies tested to date can be found at the HMS LINCS website (http://lincs.hms.harvard.edu/resources/). The extent of validation is quantified on a level between 0 and 2: "Level 0" represents antibodies for which staining was not detected using tissues for which the antigen is thought to be present based on published data; "Level 1" represents the expected pattern of positive staining in a limited number of tissues types (e.g. CD4 antibody in tonsil tissue alone); "Level 2" represents the expected pattern of positive staining in all tissues or tumor types been tested (N>= 3). Higher levels will be assigned in the future to antibodies that have undergone extensive validation; for example, side-by side comparison of against an established IHC positive control.

Mounting & de-coverslipping
Immediately prior to imaging slides were mounted with 1x PBS or, if imaging was expected to take longer than 30 minutes (this occurs in the case of samples larger than 4 cm 2 corresponding to about 200 fields of view with a 10X objective) with 1x PBS containing 10% Glycerol. Slides were covered using 24 x 60mm No. 1 coverslips (VWR 48393-106) to prevent evaporation and to facilitate susbequent decoverslipping via gravity. Following image acquisition, slides were placed in a vertical staining jar containing 1x PBS for at least 15 min. Coverslips were released from slides (and the tissue sample) via gravity as the slides were slowly drawn out of the staining jar.

Fluorophore inactivation
After imaging, fluorophores were inactivated by placing slides horizontally in 4.5% H 2 O 2 and 24 mM

Image processing
Quantitative analysis of tissue images is a complex, in large part because cells are close packed.
Background can be uneven across large images and signal to noise ratios relatively low, particularly in the case of anti-phospho-protein antibodies. We have only started to tackle these issues in the case of high dimensional t-CyCIF data and users of the method are encouraged to thoroughly research image processing methods themselves. We expect to update methods and algorithms described here on a regular basis and users are encouraged to check our Web site for additional information.

Background subtraction and image registration
Background subtraction was performed using the previously established rolling ball algorithm (with a 50-pixel radius) in ImageJ. Adjacent background-subtracted images from the same sample were then registered to each using an ImageJ script as described previously 1 . In brief, DAPI images from each cycle were used to generate reference coordinates by Rigid-body transformation. To generate virtual hyper-stacked images, the transformed coordinates were applied to images from four channel imaging of each CyCIF cycle.

Single-cell segmentation & quantification
To obtain intensity values for single cells, images were segmented using a previously described 2 conventional waterfall algorithm based on nuclear staining by Hoechst 33342. Images were binarized in the Hoechst channel and then converted into regions of interest (ROIs) for each cell. The Watershed algorithm in ImageJ was then applied to enlarge ROIs (by 3 pixels in the case of 10 x images) and encompass a significant portion of the cytoplasm and membrane for each cell. ROIs were then used to compute intensity values from all channels. All scripts can be found in our Github repository (https://github.com/sorgerlab/cycif).

High-dimensional single-cell analysis by t-SNE and EMGM
Raw intensity data generated from registered and segmented images was imported into Matlab and converted to comma separated value (csv) files. The viSNE implementation of t-SNE and EMGM algorithms from the CYT single-cell analysis package were obtained from the Pe'er laboratory at Columbia University 3 . Intensity-based measurements (such as flow cytometry or imaging cytometry) of protein expression have approximately log-normal distribution 4 , hence, t-CyCIF raw intensity values were first transformed in log or in inverse hyperbolic sine (asinh) using the default Matlab function or the CYT package 3 , respectively. Between-sample variation was then normalized on a per-channel basis by aligning intensity measurements that encompass values between 1 st and the 99th percentile (using the CYT package). Data files were aggregated and used to generate viSNE plots. All viSNE/t-SNE analyses used the following settings: perplexity -30, epsilon =500, lie factor =4 for initial 100 iterations and lie factor -1 for remaining iterations. For EMGM clustering, k = 8 was used for dividing samples into groups. Intensity values from all antibody channels (plus area and Hoechst intensity) were used for unsupervised clustering.

Calculating Shannon entropy values
Images were divided into regular grids and 1000 cells from each region used to calculate the nonparametric Shannon entropy as follows: where E1(s) is the Shannon entropy of signal s; s i is the per-pixel intensity of signal s at a given point. For GBM composites images ( Figure 5, Figure 6 & Table S5) https://omero.hms.harvard.edu/webclient/?show=dataset-2041      Table S2. List of antibodies used to stain PDAC sample presented in Figure 3 and TMA in Figure 4.