Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Deep Profiling of Mouse Splenic Architecture with CODEX Multiplexed Imaging

Deep Profiling of Mouse Splenic Architecture with CODEX Multiplexed Imaging Resource Deep Profiling of Mouse Splenic Architecture with CODEX Multiplexed Imaging Graphical Abstract Authors Yury Goltsev, Nikolay Samusik, Julia Kennedy-Darling, ..., Gustavo Vazquez, Sarah Black, Garry P. Nolan Correspondence gnolan@stanford.edu In Brief A DNA barcoding-based imaging technique uses multiplexed tissue antigen staining to enable the characterization of cell types and dynamics in a model of autoimmune disease. Highlights d Autoimmunity analyzed by multiplexed DNA-tagged antibody staining (CODEX) d CODEX data reveal pairwise interactions and niches changing with disease d First tier of neighbors significantly impacts marker expression in the index cells d Changes in splenic morphology correlate with shifts in cell frequencies Goltsev et al., 2018, Cell 174, 968–981 August 9, 2018 ª 2018 The Author(s). Published by Elsevier Inc. https://doi.org/10.1016/j.cell.2018.07.010 Resource Deep Profiling of Mouse Splenic Architecture with CODEX Multiplexed Imaging 1,2 1,2 1 1 1,3 1 Yury Goltsev, Nikolay Samusik, Julia Kennedy-Darling, Salil Bhate, Matthew Hale, Gustavo Vazquez, 1 1,4, Sarah Black, and Garry P. Nolan * Department of Microbiology and Immunology, Baxter Laboratory in Stem Cell Biology, Stanford University School of Medicine, Stanford, CA, USA These authors contributed equally Present address: Smart Tube, Inc., 922 Industrial Avenue, Palo Alto, CA 94303, USA Lead Contact *Correspondence: gnolan@stanford.edu https://doi.org/10.1016/j.cell.2018.07.010 SUMMARY tation and emission spectra makes it hard to image more than 4–5 fluorophores with conventional fluorescent microscopy. A highly multiplexed cytometric imaging approach, Yet considerably more surface markers are needed for precise identification of cellular subsets and their activation state (Chat- termed co-detection by indexing (CODEX), is used topadhyay and Roederer, 2012). Approaches have been devel- here to create multiplexed datasets of normal and oped to overcome such limitations (Schubert et al., 2006; Gerdes lupus (MRL/lpr) murine spleens. CODEX iteratively et al., 2013), but these protocols have required multiple stain/ visualizes antibody binding events using DNA barc- strip/wash cycles of the antibodies that can be time consuming odes, fluorescent dNTP analogs, and an in situ or lead to sample degradation over the iterations. polymerization-based indexing procedure. An algo- The technique described here (CODEX, for CO-Detection by rithmic pipeline for single-cell antigen quantification indEXing) extends deep phenotyping capabilities of flow and in tightly packed tissues was developed and used mass cytometry (Spitzer et al., 2015; Bendall et al., 2011)to to overlay well-known morphological features with most standard three-color fluorescence microscope platforms de novo characterization of lymphoid tissue architec- for imaging of solid tissues. Accurate highly multiplexed single- ture at a single-cell and cellular neighborhood levels. cell quantification of membrane protein expression in densely packed lymphoid tissue images (which was once deemed We observed an unexpected, profound impact of the impossible [Gerner et al., 2012]) was achieved using polymer- cellular neighborhood on the expression of protein ase-driven incorporation of dye-labeled nucleotides into the receptors on immune cells. By comparing normal DNA tag of oligonucleotide-conjugated antibodies, combined murine spleen to spleens from animals with sys- with an image-based expression estimation algorithm. Auto- temic autoimmune disease (MRL/lpr), extensive matic delineation of cell types from multidimensional marker and previously uncharacterized splenic cell-interac- expression and positional data generated by CODEX enabled tion dynamics in the healthy versus diseased state deep characterization of cellular niches and their dynamics dur- was observed. The fidelity of multiplexed spatial cy- ing autoimmune disease both for major and rare cell types popu- tometry demonstrated here allows for quantitative lating mouse spleen. A rich source of multivariate data has been systemic characterization of tissue architecture in generated and provided for the community to further efforts in normal and clinically aberrant samples. developing approaches for image analysis, tissue architecture mapping, and rare cell-type detection. INTRODUCTION RESULTS Dramatic immune tissue re-organization has been seen in lupus erythematosus, where a variety of organs (from skin, to kidney, Single Base Primer Extension Enables Multiplexed and other body organs) can be targeted in relapsing-remitting Antigen Staining flares. One example of such reorganization is pronounced DNA provides an ideal substrate for designing molecular tags lymphadenopathy and splenomegaly observed in lupus models due to its combinatorial polymer nature. An indexable tagging (Lieberum and Hartmann, 1988; Jacobson et al., 1995). Using system whereby tags are iteratively revealed in situ by a stepwise mice with MRL/lpr genotype (Kanauchi et al., 1991), we sought enumeration procedure was designed. Antibodies (or other to systematically characterize microenvironment and cell inter- affinity-based probes) are labeled with uniquely designed oligo- actions associated with changes in immune organ architecture nucleotide duplexes with 5 overhangs that enable iterative step- and the progression of autoimmune disease. To this end, we wise visualization (Figure 1A; Video S1, part 1). Cells are stained devised a multiplexed microscopy technique that allows a pre- with a mixture of all tagged antibodies at once. At each rendering cise mapping of cell types in tissues. Significant overlap in exci- cycle, the cells are exposed to a nucleotide mix that contains one 968 Cell 174, 968–981, August 9, 2018 ª 2018 The Author(s). Published by Elsevier Inc. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Figure 1. Sequential Primer Extension on Samples Stained with DNA-Barcoded Antibodies Enables Unlimited Level of Multiplexing (A) CODEX schematic diagram. (B and C) Mouse spleen cells were fixed and co-stained with conventional TCR-b Ax488 antibody and CD4 antibody conjugated to CODEX oligonucleotide duplex as in first round of (A). After staining, cells were incubated in extension buffer with dG and dUTP-Cy5 either without (B) or with (C) Klenow exo- polymerase. Note that TCR-b-positive T cells in (B) and (C) are indicated by Ax-488 staining. Dependent upon the addition of Klenow, TCR-b-positive CD4 positive T cells are seen as a Cy5 positive subset of TCR-b-positive T cells in (C). (D) Spleen cryosection stained with B cell-specific B220-APC (red) and T cell-specific TCR-fluorescein isothiocyanate (FITC) (green) show mutually exclusive staining pattern in the marginal area between B cell follicle and the white pulp. (E) Spleen cryosection stained with CODEX DNA-tagged B220 (red) and CODEX DNA-tagged TCR-b (green) shows staining similar to the one observed with regular antibodies in (D). (F) Spleen sections were co-stained with regular B220-FITC and two antibodies (ERTR7 and CD169) tagged with cycle 1 CODEX DNA duplexes. Localization of marginal zone CD169 positive macrophages in the area between the ERTR7 positive splenic conduit of the white pulp and the B220 positive follicular B cells (D) as reported previously has been observed. See also Figure S1 and Video S1, part 1. of two non-fluorescent ‘‘index’’ nucleotides and two fluorescent is paused at the indexing position by omitting one of the indexing labeling nucleotides. The index nucleotides fills in the first index (walking) bases from the labeling mix (as done in this study) or position across all antibodies bound to the cells. However, the potentially by use of reversible terminators (Video S1, part 2). DNA tags are designed such that only the first two antibodies Importantly, the system enables multiplexed tissue imaging are capable of being labeled with one of the two fluorescent analysis by means of a standard fluorescence microscope. dNTPs—and only if the index nucleotide was previously incorpo- To test the premise of the system, isolated mouse spleen cells rated. Those two antibodies are then imaged by standard fluo- were incubated with a CD4 antibody conjugated to an indexing rescence microscopy. Then the fluorophores are cleaved and oligonucleotide duplex (as represented by Ab1 in Figure 1A). In washed away, and the sample is ready for the next cycle where this trial experiment, T cell receptor b (TCR-b)-Alexa 488 was a different indexing nucleotide is used. At the end of the multi- used as a counterstain. A single round of primer extension was cycle rendering protocol each pair of antibodies is visualized at done with a mix of unlabeled dGTP and dUTP-ss-Cy5. A cell a known, pre-defined cycle of the indexing protocol, and the population positive for both CD4 and TCR-b was observed by multiparameter image can be reconstructed. The polymerase flow cytometry. Observation of this population was dependent Cell 174, 968–981, August 9, 2018 969 on the addition of Klenow DNA polymerase to the reaction (Figure 2). Indeed, use of this compensation method resulted mixture (Figures 1B and 1C) indicating specific antibody in a considerable (approximately 2-fold) reduction of spillover rendering by primer extension. Similarly, in tissue sections, signal (especially pronounced for CD4/CD8a co-distribution— CODEX tag-conjugated antibodies produced lineage-specific Figures 1L and 2E). staining comparable to regular fluorescent antibodies (compare A 30-antibody panel was therefore designed to identify staining patterns of B220-CODEX and B220-APC in mouse splenic-resident cell types (lymphocytes, macrophages, micro- spleen, Figures 1D–1F). vessels, conduit system, and splenic stroma; Figure 3A; Table A simulated multicellular mix was produced by combining 30 S1) and applied to the cryosections of spleens from wild-type aliquots of mouse splenocytes barcoded with pan-leukocytic (3 spleens) and MRL/lpr mice (6 spleens) (Figure S3). Four major CD45 antibody labeled with one of 30 distinct CODEX tags (Fig- classic splenic compartments: red pulp, B cell follicle, PALS ure S1A). The visualization of the CODEX 15-cycle staining (periarteriolar lymphoid sheath), and marginal zone (MZ) (Fig- pattern showed even cycle-specific signals, with low back- ure 3B) could be easily discerned in CODEX imaging data (Fig- ground (average signal to noise 85:1), efficient (98%) release ure 3A). A total of 734,101 segmented cells were identified, of fluorophore by inter-cycle TCEP (Tris(2-carboxyethyl)phos- and by means of X-shift clustering (see STAR Methods) their phine hydrochloride) cleavage and no signal carryover between expression profiles were grouped into 27 broadly defined pheno- cycles (Figures S1B–S1D and S1F). Linear regression analysis typic groups (Figures 3C–3F; Table S2) most of them matching revealed low signal deterioration (at 0.79% per cycle) and to known cell types. Compared to CyTOF data on splenocytes acceptable background (starting at 1.1% and increasing at isolated from non-enzymatically homogenized spleen, CODEX 0.06% per cycle, Figures S1E and S1F). in situ analysis produced a similar distribution of cell counts for To further reduce the signal loss associated with accumulation major cell type. Yet being a non-disruptive technique CODEX of polymerization errors and to allow larger panels without identified larger numbers of resident and stromal cell types increasing the length of tagging oligonucleotides, an approach such as erythroblasts and F4/80 macrophages than CyTOF did based on primer-dependent subpanels was devised (Fig- (Figure 3E). Notably even rare computationally derived cellular hi – hi + ure S2A). The feasibility of this design expansion was tested by phenotypes (e.g., CD4 /CD3 /MHCII cells and CD11c B cells) staining mouse splenocytes with a 22-plex set of antibodies. closely matched the cell types previously observed in murine Each of the antibodies was conjugated to the three versions of spleen (LTi cells [Robinette et al., 2015] Figures S4A–S4C, the CODEX duplex tag—with same terminated top oligonucleo- S5C, S5F, and S5I and age-associated B cells [ABCs], Figures tide and three kinds of the tagging oligo (Figure S2B). Thus, every S4D, S4E, S5B, S5E, and S5H). antigen was detected thrice (bringing the overall number of detections to 66) and only after annealing of a panel-specific Pairwise and Combinatorial (i-niches) Statistics of Cell- activator oligonucleotide. We found that the signal for same to-Cell Contacts in Mouse Spleen antibody was consistent across the three primer-dependent To provide a high-level view of the cell-type interaction land- batches (Figure S2C). Thus, panel-activator design extends scape, the total counts of contacts between every pair of cell CODEX to a theoretically unlimited multiplexing capacity, types in the Delaunay neighborhood graph (Gabriel and Sokal, bounded only by the speed and resolution of the imaging pro- 1969)(Figure 4A and associated Mendeley dataset) for each cess itself and the time required for each imaging cycle. condition was determined. The specificity of cell-to-cell interac- tion was estimated from the ‘‘log odds ratio’’ metric (log-ratio of Benchmarking CODEX observed and expected probabilities of contacts between 2 cell To validate the quantitative performance of CODEX, cells freshly types) (Table S3). When visualized as heatmaps, this metric re- isolated from mouse spleens were co-analyzed by mass-cytom- vealed a significant non-random distribution of cells in the etry (CyTOF) and CODEX using identical 24-antibody panels spleen. In the majority of cases, cell types were either selectively (Table S1). Use of the same antibody clones and the same sple- associating or avoiding each other (red or blue on the heatmap) nocyte preparation ensured the validity of comparisons. Antigen pointing to prevalence of specific cell-to-cell interactions in co-expression signals from CODEX, obtained from image seg- shaping the spleen architecture. The major splenic anatomic mentation (see STAR Methods), were consistently similar to compartments were reflected in two large mutually exclusive CyTOF (Figure 2B, for direct comparability, both CODEX and clusters of positive associations, which appeared to correspond CyTOF data are plotted on a linear scale). to red pulp and the white pulp, respectively (indicated with black Consecutively, CODEX was applied to tissue sections. In rectangular outlines on Figure 3G). For example, a significant contrast to dissociated cells spreads (Figure 2A), cells in tissue positive association was observed between F4/80 macro- sections are adjacent to each other—with large fractions of phages and erythroid cells, as these cell types are both found membranes in direct contact (Figure 2C). Therefore, neighboring in the red pulp and are closely associated in so-called erythro- cells can contaminate each other’s signals during the quantifica- blast islands (Socolovsky et al., 2007). An avoidance of interac- tion phase (Figure 2D). To address this latter challenge, a novel tion was observed between T and B cells, reflecting concentra- linear algorithm for 3D positional spillover compensation was tion of these cell types in B cell follicles and PALS, respectively created. This algorithm is based on the same principles used (Figure 3G). in fluorescent spillover compensation in traditional flow cytome- Unexpectedly, a consistently high degree of association was try, except that our algorithm performs compensation between observed between the cells of the same phenotypic class physically adjacent cell based on their surface contact ratios (Figure 3G, red diagonal), suggesting that homotypic adhesion 970 Cell 174, 968–981, August 9, 2018 Figure 2. Accuracy of Surface-Marker Quantitation by CODEX (A) Microscopic image of mouse splenocytes stained with a 24-color antibody panel, showing one cycle of CODEX antibody rendering. Cell contours show the outlines produced by the cell segmentation algorithm (B) Comparison of single-cell expression data derived from dissociated mouse splenocytes on an identical 24-color panel using CODEX and CyTOF. (C) Example segmentation in a mouse spleen section based on combining nuclear and membrane (CD45) channel. (D) Graphical explanation of the algorithm for compensating the spillover between neighboring cells using a cell-by-cell compensation matrix. (E) Biaxial plots of segmented CODEX data acquired in mouse (BALBc) spleen sections. The presence of double-positive cells in the upper quadrant is used as an estimate of lateral signal bleeding explained schematically in (D). Three combinations of mutually exclusive lineage markers are shown to demonstrate the range of effect of the compensation algorithm on reduction of lateral signal bleed. See also Figures S4 and S5. constitutes a major force driving the architecture of immune tis- previously. We define here an indexed ‘‘niche’’ (i-niche) as a ring sue. This observation held true both for the major constituents of of cells (excluding the central, or here defined as ‘‘index’’ cell) in white pulp, T and B cells, as well as for rare cell types such as no specific circumferential order that are Delaunay neighbors of natural killer (NK) cells. Interestingly, even though CD8 and the index cell (Figure 4A and STAR Methods). We identified 100 CD4 T cells tended to mix in the PALS, their mutual distribution of the major i-niches (by K-means clustering) according to the was nonrandom and consisted of intertwined threads of homo- relative frequency of the identified cell types present in the ring typic cells (Figure S4H). Interestingly, as an aside, similar struc- of cells surrounding the index cell (Figures 4A, 4B, and 4F). tures could be reproduced in vitro by incubating heterotypic Most i-niches could be readily mapped into one of major mixtures of sorted splenic cell populations (Figures S4I and anatomical compartments of the spleen (B cell follicle, PALS, S4J). These data suggest that homotypic cell association marginal zone, or red pulp, per Figure 4G). In most cases, any might be an important driver of the white pulp substructure given i-niche resided within a single anatomical compartment and is worth further investigation. (although several i-niches were observed in more than one The precision in situ cytometry analysis of CODEX data al- compartment), and every splenic compartment was populated lowed enumeration of cellular contexts in a manner not possible by many i-niches (Figure 4G). Cell 174, 968–981, August 9, 2018 971 Figure 3. CODEX Analysis of Mouse Spleen Cryosections Co-stained for 28 Antigens (A) Three collated images on the left correspond to the legend of antibody renderings per cycle, gross morphology photograph of MRL/lpr (left), and normal (right) spleen embedded in optimal cutting temperature (OCT) block prior to sectioning. Green corresponds to antibodies rendered by extension with dUTP-Cy5; red, dCTP-Cy3 On the right, collage of the CODEX multicycle data for normal spleen (BALBc-2) and early MRL/lpr spleen (MRL/lpr 4). All images are derived from a single scan with a 403 oil objective of an area covered by 63 tiled fields. (B) Schematic diagram of major known splenic anatomical subdivisions drawn based on cell distribution in BALBc-1 replicate. (C) An exemplary profile of Vortex cluster (B cells) used for manual matching of clusters to known cell types. (D) Minimal spanning tree (MSP) built for all clusters identified by Vortex analysis. On the left middle and right panels, the MSPs are colored by expression levels of B220, TCR, and CD71 accordingly to indicate location of B cells T cells and erythroblasts on the tree. (E) Circle chart showing for several major cell types their fraction of total cells as identified by CODEX analysis of splenic tissue and CYTOF analysis of isolated BALBc splenocytes (F) Post-segmentation-derived diagram of identified objects (cells) colored according to cell types in BALBc-1 replicate. Full-size diagrams for every tissue analyzed in this study are available online (see STAR Methods) (G) Average cell-type to cell-type interaction strength heatmap for BALBc samples. Color from blue (<0) to white (around 0) to red (>0) indicates log of odds ratio of interaction (ratio of observed frequency versus expected frequency of interaction). The rows and columns are in the same order (annotation on the right). Black outlines indicate two largely exclusive mega-clusters of cross-interacting cell types loosely matching the cell types populating the red and the white pulp. See also Figure S5 and Video S2. 972 Cell 174, 968–981, August 9, 2018 Figure 4. Unbiased Identification of i-niches in Multidimensional CODEX Data (A) On the left diagram explaining the terminology used for defining i-niche (a ring of first tier neighbors for central cell). On the right Delaunay triangulation graph used for identification of first tier of neighbors for every cell. (B) Heatmap depicting frequency of cell types in 100 types of i-niches identified by K-means (K = 100) clustering of all index cells in the dataset (each cell is an index cell for its i-niche) based on frequency of different cell types in the first tier of neighbors. The color indicates the average fraction of corresponding cell type in the the i-niche. (C) An example of marginal zone and follicular (B-zone) B cells defined by residence in distinct i-niches (e.g., marginal zone i-niche includes a marginal zone macrophage marked by letter H and green color). Positions of B cells in each i-niche is marked with red circles over the schematic of BALBc spleen. (D and E) Two heatmaps from top to bottom show average expression of selected surface markers measured in a central cell across 100 i-niches (same left to right order as in B) when central cell is B cells (D) or CD4 T cell (E) accordingly. The color indicates the relative level of surface-marker expression as measured across dataset. Gray columns indicate absence of cells in corresponding niches. Two orange rectangles over top heatmap indicates position of i-niches with high CD35 (containing FDCs and marginal zone macrophages). Cyan rectangle shows location of family of i-niches with high content of F4/80 macrophages and low B220 and CD19 in central B cells. Purple rectangle indicates family of i-niches enriched with ERTR-7 positive stroma. Below top heatmap, location of selected i-niches shown in (E) are indicated. Over bottom heatmap, yellow rectangle indicates the family of i-niches with dominating presence of B cells. Two green rectangles indicate family of niches with high levels of CD90 and CD27 in the index CD4 T cells. (F and G) Abundance of 100 i-niches in normal spleen (top bar graph) (F) and relative distribution of i-niches (G) between splenic histological subdivisions (PALS, red pulp, marginal zone, and B-zone) shown as a heatmap. To illustrate a variety of tissue distribution pattern by i-niches an overlay of selected i-niches over a schematic of normal spleen (BALBc-1) is shown. Heatmap color indicates fraction of corresponding i-niche per splenic anatomic subdivision. (H) Top right shows a biaxial plot of flow data for CD79b and B220 measured in isolated splenocytes. Top left shows levels of CD79b and B220 in central B cells as measured across all 100 i-niches. To illustrate i-niche-dependent variability of surface-marker expression, images of central cells (marked with red cross) with levels of surface marker indicated in pseudocolor palette are shown for selected exemplary i-niches in the bottom panels. See also Figures S4K and S4L. Cell 174, 968–981, August 9, 2018 973 In our definition, the index cell in the center can be any cell. levels on splenic B cells could be, to a large degree, accounted Thus, i-niches enabled subsetting the common cell types based for by the niche composition around those B cells—and that on cellular context (microenvironment). For example, B cells sur- the expression levels on these cells might be influenced by (or rounded by only B cells (red arrow Figure 4B, i-niche #96) can be influences) the cells in their immediate surrounding. seen primarily in the follicular zone B cell region (Figure 4C, left The overall utility of the i-niche in determining any given sur- panel), while presence of CD169 positive marginal zone macro- face-marker expression value for an index cell was evaluated phages mapped the B cells in such i-niches to the marginal zone by constructing a linear regression model of marker expression (red arrow Figure 4B, i-niche #33, Figure 4C, right panel). In the using both the cell-type identity and the i-niche constituency in case of T cells, CODEX data enabled precise selection of an a two-featured variable model (the other variable being the cell- important migratory subset of T cells known to be residing in type identity). Notably, adding the i-niche information as a depen- ERTR7 enriched niches (Burrell et al., 2015) (in Figures 4B and dent variable significantly improved the fitness of the model for all 4E see the purple rectangle indicating a family of niches where markers (Table S4) with highest improvement F-values for CD90, index T cells contact ERTR7 stroma; as well as Figures S6A B220, CD21/35, and ERTR7 and the lowest prediction rates for and S6B). Taken together, we see that surface-marker expres- Ly6G, CD5, CD11b, CD5, and TCR-b. Thus, the high variability sion alone is insufficient to associate many cell subsets with a in B220 expression levels are highly related to the i-niche in which given tissue subcompartment (e.g., CD4 T cells can be found the B cell resides. In other words, B220 expression levels can be both in the PALS and in the red pulp). However, i-niche designa- location specific and are dependent on the i-niche partners. As a tion does provide such mapping data (most of i-niches were en- counter point, the data also show that i-niche does not reliably riched within a specific splenic subdivision Figure 4G) and as predict expression of other proteins, such as CD5 or TCR-b, such T cells associated with ERTR7 positive stroma in fact the expression levels of these receptors is relatively constant localize primarily to the red pulp (Figures S6A and S6B). This rai- across the i-niches (Figure 4E). This result quantitatively demon- ses interesting questions—can new cell types, or functional sub- strates that the i-niche (neighbors) determines a significant pro- sets, be discerned by this approach? What is the frequency of a portion of variance in the expression of certain markers. Overall, repeated i-niche structure that must be observed to suggest a we observed that that many splenic cell types populate a wide va- function? And what would constitute a proof that a given i-niche riety of i-niches (Figures S4K and S4L), suggestive of a multiplicity corresponds to a new cell type or functionality? of functional state for any given immune cell type. Further, tissue locale (i-niches) is a powerful indicator of potential differential Cell-Surface-Marker Expression Depends on Local function (to the extent tissue locale drives function), and these Neighborhood deterministic changes in surface marker protein expression are One approach to address these latter questions is to consider surrogate indicators of this locale or function. the phenotypes of the index cell in various i-niches. We observed that for several index cell types (e.g., for B and T cells) there was Changes in Splenic Composition Associated with significant biasing of the surface-marker expression depending Disease Progression on the i-niche in which the index cell resides (e.g., see selected A comparable region of spleen was visualized by CODEX for 3 cases marked with colored rectangles above the heatmaps in normal BALBc spleens, and 6 spleens from MRL/lpr mice. Image Figures 4D and 4E). To assure that it was not a quantitation arti- segmentation revealed strong variation in cell counts between fact, we mapped back to the image the index B cells where the the norm and the disease (Figure 5B) for most (19 out of 27) of levels of CD79b (a co-activator chain of the B cell receptor com- the cell types identified by X-shift clustering. Examples include plex) and B220 (a splice isoform of CD45 membrane phospha- a dramatic increase in CD71 erythroblasts (green cells on Fig- tase) would be niche dependent. We found that the index B cells ure 5A maps), a reduction in numbers of B cells and FDCs (follic- int lo + that were B220 , CD79b (i-niche ‘‘59’’) resided on the ular dendritic cells), and increases in so-called B220 DN T cells boundary between the PALS and the follicles (Figure 4H, image (CD4/CD8 double-negative B220 T cells), which have been pre- lo montage on the bottom). Index B cells that were B220 , viously characterized as a hallmark of the MRL/lpr progression int CD79b (i-niche ‘‘91’’) were mostly found in the red pulp. And, (Koh et al., 1995) which could also be identified by fluores- int/hi hi B cells that were B220 , CD79b (i-niche ‘‘76’’) were yet cence-activated cell sorting (FACS) (Figures S1F and S1G), different again and were found at the boundary of the red pulp thus ruling out the possibility that this unusual cell type being a and the follicles. When measured in cell suspensions, CD79b result of image segmentation errors. These and other changes is co-expressed with B220 at various ratios (see CyTOF plot sin- were used to broadly classify the MRL/lpr spleens into early, in- gle-cell splenocytes Figure 4H, top-right panel). Such a distribu- termediate, and late disease stages (Figure S3). tion of expression is sometimes attributed to staining variability, Despite consistent presence of ‘‘homotypic interactions’’ di- measurement noise, or a simple lack of understanding of the agonal and larger cell-adjacency clusters corresponding to red underlying biology. However, as seen on Figure 4H, upper-left and white pulp in odds ratio heatmaps across disease (Fig- panel, there is a non-random pattern of CD79b and B220 expres- ure S3), a deeper statistical analysis revealed many disease- sion across the central cell of the corresponding i-niches, and, associated changes in frequency of contacts between cell types depending on the B220/CD79b levels, the i-niches (the central (see Table S3). Among the changes we observed an increase in – + cells) map to specific regions in the splenic architecture (Fig- interaction between B cells and CD4 /CD8 conventional den- ure 4H, lower panel). These observations suggest that the dritic cells (cDCs) in the early MRL/lpr spleen compared to spread of the CD79b-B220 levels as well as of other marker normal (Figure 5C, left panel), suggesting an increase in B cell 974 Cell 174, 968–981, August 9, 2018 Figure 5. Autoimmune Disease Drives Changes in Splenic Composition and Cell-to-Cell Interactions (A) Post-segmentation diagrams of all objects (cells) colored according to cell types (see color map in Figure 3F) for all normal and MRL/lpr tissue sections imaged in the study. Full-size diagrams are available for every tissue analyzed in this study are available online (see STAR Methods). (B) Stacked bar graphs show dynamics of cell counts across dataset for manually annotated Vortex clusters (cell types on the left) across progression from normal to afflicted spleen. Colored bar sections indicate fraction of the total cells as detected at a particular stage/samples (1–9 annotation on the top). Cell types were split into four types according to the dynamics of counts across dataset as represented by average relative (normalized to 1) count; see line graphs on the right; x axis corresponds to stage/sample id. (C) Two examples of change in cell-to-cell interaction frequency during disease progression between the B cells and dendritic cells in normal and early MRL/lpr spleen and between B220 DN T cells and CD4 T cells during progression from early MRL/lpr to intermediate. (D) Co-distribution of odds ratio log fold [log(odds ratio in early MRL/lpr) – log(odds ratio in BALBc)] on x axis and change in counts of interactions for early MRL/lpr versus control (BALBc) comparisons (on y axis). (E) Co-distribution of cumulative cell-frequency change [celltype1 freq. change + celltype2 freq. change] on x axis and change in counts of interactions for early MRL/lpr versus control (BALBc) comparisons (on y axis). (F) Bar graph showing chi-square values across conditions computed for odds ratio and direct interaction counts. See also Video S2 and Figure S6. activation. We also observed a higher interaction frequency of S6G). In the intermediate- and late-stage MRL/lpr spleens, there granulocytes with T cells (Figure S6C), dendritic cells (Fig- was a significant increase in interaction of B220 DN T cells with + + ure S6D), and erythroblasts (Figure S6E), and a higher number CD4 T cells (Figure 5C, right panel), CD8 T cells, erythroblasts, of contacts between erythroblasts and various kinds of stromal and a variety of other cell types compared with numbers of these cells, as well as B220 DN T cells (Table S3; Figures S6F and interactions in the early MRL/lpr stage (Table S3 and Figures Cell 174, 968–981, August 9, 2018 975 S6G–S6I, see online resource, STAR Methods, for more exam- For the rest 24 of the 26 changing interactions mentioned ples). So, while there was no obvious gross rearrangement of above at least one of the cells of the pair was scored as signifi- the tissues, many homotypic and heterotypic cell-cell associa- cantly (FDR <0.05) changing the frequency across scored condi- tions were altered, prompting a key question: what are the tions (Table S3, last column of the ‘‘EarlyMRL vs BALBc control’’ main factors driving this disruption? spreadsheet). We therefore conclude that—at least in the diseased state of early-stage MRL/lpr—most of the change in Disease-Driven Change in Cell Counts Determines the counts of cell-cell interactions are driven simply by increases Frequency of Specific Cell-to-Cell Contacts or decreases in cell-type frequencies. What could be the drivers of changes in frequency of pairwise As an additional evidence, c statistics were used to compare cell-cell contacts? If the kinetics of pairwise cell contacts follows the total magnitude of changes in pairwise cell-type interaction a rate law, one possibility would be that modulation of specific matrices (total interaction count) versus changes in log-odds cell-to-cell interaction potential—or ‘‘attraction’’ (for which the ratio matrices (propensity for non-random interaction). The c odds ratio score was used as an estimate across this study)— deviation (sum of squares of Z-score-normalized values) was is the main driver. In other words, it would be expected that, computed for each disease matrix compared to the control. when the affinity of such an interaction goes up, the fraction In every case, the c values of cell-interaction matrices were of interacting cells of a given cell pair would increase. At larger than of the respective log odds ratio matrices of the the same time, even in the absence of change in cell-to-cell same biological sample (Figure 5F). This suggests that as affinity, the absolute number of the cell-cell pairs (defined here the cell-type frequencies change due to disease progression, the absolute numbers of interactions change dramatically as cell pair aggregates, or CPAs) and the number of interacting cell pairs should correlate with the frequencies of interacting whereas the frequency-normalized likelihoods of cell interac- cell types (analogous to concentrations in the rate law equation). tions change to a much smaller extent indicating a great degree Importantly, the latter scenario could be as biologically signifi- of robustness of the ‘‘design principles’’ of the splenic tissue and cant as the former. Finally, some of the cell-type contacts may that many of the more dramatic disease-associated variations be observed due to low cellular motility of randomly meeting occur primarily through the shift in cell numbers. cells. Such interactions would not produce spatially defined This analysis implies that the degeneration of the tissue integ- sub-splenic CPAs and would have and odds ratio close to 1. rity in the MRL disease largely follows dramatic changes in cell- The perturbation introduced to normal splenic composition type frequencies. At the same time, there were notable excep- with MRL/lpr genotype allowed us to examine the mechanisms tions from this trend, where the changes in observed cell-type implicated in transition from normal to diseased spleen. In short, pairing frequencies could be largely explained by shifts in the we found that, for most cell-cell pairs observed, the mutual cell-type interaction likelihoods (log-odds ratios). While further attraction (as quantified by the interaction log odds ratio) was work is required to determine which of these changes are instru- not the primary determinant driving the change in counts of inter- mental to the MRL disease state, our findings suggest that such acting cell pairs between the MRL/lpr and the norm. In Figure 5D, differential analysis can be applicable in other diseases and, we plot the change in counts of interactions of two cell types possibly, could be used to discover cell-type interactions that (e.g., A:B) between the MRL/lpr and the normal BALBc spleens. are targetable from a therapeutic standpoint. Each dot represents a pair of cell types. The value on the y axis is the difference in the total number of observed interactions be- Reorganization of Cells in Disease-Associated Tissue tween BALBc and MRL/lpr. The x axis shows the difference be- Substructures tween log odds ratios of interactions between the same condi- We cataloged the cell-cell interaction ‘‘connectivity’’ in a circular tions. There was no overall correlation observed (R = 0.058). correlation diagram. Rarely, if ever, there was any cell type found In contrast, we observed a correlation (R = 0.288) between adjacent to only one other type of cell. The highest degree of the cell count changes and the interaction changes (Figure 5E). connectivity was observed for the most abundant cell types In agreement with those observations, we saw that out of the such as B cells in normal spleen and erythroblasts (Figure 6A) 26 top-scoring (false discovery rate [FDR] <0.05 and change in early MRL/lpr. This high connectivity in turn led to large effect in absolute interaction counts >150) cell-type pairs of this on i-niches caused by changes in cell numbers associated with cross comparison, only 2 showed corresponding significant progression of disease from normal to autoimmunity. The most (FDR <0.05) change in odds ratio score. Curiously, these two in- dramatic changes in cell frequencies were the increase in eryth- teractions with a modest 1.5 times increase in interaction count roblasts in the early MRL/lpr and the emergence of B220 DN and, concomitantly, a 0.8 increase in log odds ratio score were T cells in late MRL/lpr—which were associated with the appear- the ones between the CD4 or CD8 T cells and ERTR7 stroma ance of novel i-niches relative to the normal spleen (spatial (see Table S3, rows 6 and 7, and Figures S6A and S6B). Visually localizations of B220 DN T cell-dominated i-niche 18, erythro- they appeared as persistent co-clustering of T cells with ERTR7 blast-driven i-niche 29, and B cell-rich i-niche 96 are shown on stroma despite the overall drop of T cell numbers in the ‘‘early’’ Figure 6B, and their cell-type composition is shown on heatmap MRL/lpr samples. Curiously, ERTR7 positive fibers of splenic on Figure 6C). A corollary to this is the question of whether the stroma as well as ERTR7 protein itself were recently shown to presence of these cells, and new i-niches dependent on these be critically involved in T cell trafficking (Burrell et al., 2015), sug- cells, somehow changed the observable biology of the cells gesting that this increase in the spatial association could be they contact? We found some examples supporting that, reflective of the T cell activation. whereby the proximity of CD4 T cells to B220 DN T leads to 976 Cell 174, 968–981, August 9, 2018 Figure 6. Differential Effect of Disease over i-niche Presence across Dataset (A) Cell-interaction networks built for BALBc early MRL/lpr and late MRL/lpr based on the number of contacts observed between two cell types (only connections with more then 150 interactions per sample are shown on the diagrams). Thickness of connection correlates with number of contacts size of the node indicates number of cells per condition. (B) Evolution of i-niche abundance across dataset. Selected three i-niches (marked above heatmap in C depicting i-niche composition) differentially represented across dataset (changing between norm and disease) are shown. Yellow circles overlaid over blank rectangles corresponding to imaged area indicate location of i-niche. (C) Top heatmap shows frequencies of B220 DN T cells, erythroblasts, and B cells in the i-niche rings. Line above top heatmap indicates the composition of i-niches 18, 29, and 96 described in (B). Color scheme is the same middle heatmap and indicates expression of selected markers when the i-niche central cell is an erythroblast—primarily to show that CD27 is not expressed on erythroblasts in the vicinity of B220 DN T cells. Bottom heatmap indicates expression of selected (legend continued on next page) Cell 174, 968–981, August 9, 2018 977 CD4 T cell activation in spleens of MRL/lpr mice: Figure 6C The neural network highlighted the regions in each multipa- shows increased levels of CD27 expression in CD4 T cells pre- rameter spleen image that corresponded to the disease state sent in i-niches dominated by B220 DN T cells (Figure 6C, red (Figure 7B), despite having seen no images from these spleens circle). during training. To investigate the specific features learned by Other cell types noticeably changed their characteristic distri- the neural network, the cell-type compositions of the regions bution and their propensity to engage, or evade, specific cell-to- identified as diseased versus those regions identified as normal cell contacts (as estimated by odds ratio score) during disease were compared. There was significant enrichment of several cell + – + progression. For example, cells of CD106 CD16/32 Ly6C types in these regions (Figure 7C). Although some cell types en- + + CD31 phenotype were randomly distributed in the red pulp riched in diseased regions, for example, B220 DN T cells, were of normal spleens but were found to aggregate in the areas present only in the diseased tissue, the most highly enriched cell + – proximal to the marginal zone of the MRL/lpr white pulp (Fig- type (CD4 /CD8 cDCs) was present in both the disease state ures S5D, S5G, and S5J). This re-distribution correlated with and the healthy state. erythroid proliferation and reduced odds ratio score for the inter- To assess the specific contextual changes recognized by the + – + + + – action of CD106 CD16/32 Ly6C CD31 and erythroblasts in neural network, the local neighborhoods of the CD4 /CD8 cDCs lupus spleens (Table S3). that the neural network found to be enriched in MRL/lpr regions were analyzed. In these neighborhoods, we observed a signifi- + – Automatic Definition of Disease-Associated Areas in cant enrichment of other CD4 /CD8 cDCs, as well as signif- + + – – Tissue Architecture icant depletion of CD106 /CD16/32 /Ly6C /CD31 stromal cells (FDR <10 ). This suggests that the neural network had identified As noted, the analysis reveals that the development of the auto- + – immune disease in mice (as exemplified by MRL/lpr lupus) is an altered context for CD4 /CD8 cDCs (distant from stromal associated with vast rearrangement of normal spleen architec- regions) as a key descriptor for the disease. Thus, the neural ture, which is likely to cause loss of cell-cell contexts normally network approach described here enabled both automatic clas- hosting the cells crucial for proper splenic function, as well as sification of samples according to disease state and an auto- the observed emergence of novel i-niches that are not found matic identification of high-dimensional regions of interest and in the normal BALBc spleen. Additionally, certain i-niches corresponding cellular niches. were sequestered to specific anatomic compartments of the spleen, which allowed us to use such i-niches as reference DISCUSSION points to quantitatively monitor high-order morphological changes. The i-niches that in normal spleen were localized to Here the feasibility of polymerase-driven highly multiplexed visu- one distinct compartment (more than 90% of central cells reside alization of antibody binding events to dissociated single cells as within a particular splenic compartment) were used to evaluate well as tissue sections (CODEX) was demonstrated and bench- the dynamics of splenic cells associated with progression of marked. Critically, CODEX enables co-staining of all antigens autoimmune disease (Figure 7A, middle heatmap). This analysis simultaneously with the staining iteratively revealed by primer confirmed the dissipation of the marginal zone starting from extension cycles wherein no diminution of epitope signal detec- early stages of MRL/lpr and revealed a progressive distortion tion was observed. A consistent performance of CODEX in co- of PALS. Curiously, depending on whether a i-niche was based detecting up to 66 antigens was demonstrated, and the ‘‘activa- on F4/80 macrophages or primarily contained erythroblasts, the tion primer’’-based extension of the system could enable a red pulp appeared to reorganize in the diseased tissue (Fig- potentially vast expansion of CODEX multiplexing capacity. For ure 7A, right heatmap), pointing to the fact that more than the current method, fresh-frozen tissue was used yet at a cost one compartment-specific niche is required to reliably trace of testing an extensively broader collection of clones we have the fate of specific anatomic compartments. In many cases, recently succeeded in adapting the procedure to formalin fixed the definition of subsets/morphological units constituting the paraffin embedded (FFPE) archival tissue (unpublished data). tissue is subjective, yet this study employed niches that were We believe this will open the large retrospective collection of algorithmically defined. Therefore, using niches as markers of FFPE samples from clinical cohorts to multidimensional cyto- morphology can quantitatively monitor the changes of high- metric analysis. order anatomic architecture. With the help of a simple automated fluidic setup (Figure S7), To automatically isolate the specific local combinations of the CODEX platform can be performed on most three-color fluo- expression patterns characteristic of the disease state, a fully rescence microscopes with a motorized stage and thereby convolutional neural network was trained to distinguish image enable conversion of regular fluorescence microscope into a patches from normal and MRL/lpr mice. The neural network tool for multidimensional tissue rendering and cell cytometry. operated by identifying, in each training image patch, the spe- CODEX completes a 30-antibody visualization in approximately cific areas that corresponded to the disease state. 3.5 hr. Modifications to the technology that increase the markers when the i-niche central cell is a CD4 T cell. The color schemes in these three heatmaps are the same as in heatmaps in Figures 4B, 4D, and 4E. Red oval outline pinpoints i-niches with elevated CD27. Note that these i-niches as indicated by top heatmap have B220 DN T cells as a prevailing component. Lower panels show examples of central cells in i-niches marked under the lower heatmap. i-niche 50 is an example of i-niche without B220 DN T cells. Central cell does not express high CD27. i-niches 42 and 44 have high frequency of B220 DN T cells and accordingly central cells express high CD27. See also Figures S4K and S4L. 978 Cell 174, 968–981, August 9, 2018 Figure 7. i-niches and Neural Nets Provide Unbiased Way for Disease Monitoring (A) Selected i-niches (green heatmap shows i-niche composition, color scheme same as in Figure 4B) were chosen based on high (>90%) presence per single histological subdivision (blue heatmap color scheme same as in Figure 4G). Abundance of these i-niches (brown heatmap, color indicates relative abundance of corresponding i-niche as measured across full dataset) was used to judge the preservation or decay of a histological splenic subdivision corresponding to selected i-niches. (B) Red color over blue rectangle indicates regions of interest (MRL/lpr-specific regions) predicted by neural network in entire spleen images. From top left, clockwise: BALBc #3, MRL/lpr #5, MRL/lpr #7, MRL/lpr #8. (C) Cell types enriched (FDR <0.1) in MRL/lpr-specific regions (in red in B) predicted by neural network. See also STAR Methods ‘‘Neural network training.’’ measurements per cycle, reduce the cycling time, faster imaging (MRL/lpr). Much like with conventional flow cytometry, CODEX methods such as light sheet microscopy, or an increased size of discerned all major cell types commonly observed in mouse the imaging the field of view offer potential opportunities for spleen. Moreover, application of X-shift phenotype-mapping al- increasing the depth and speed of the visualization process. gorithm (Samusik et al., 2016) tailored to parsing the multidimen- Given the low cost of converting a scope to this platform sional single-cell data enabled the detection of rare cells types hi hi + (enabled by a simple fluidics device for automated sample (such as CD4 MHCII [Lti] cells, CD11c B cells [ABC cells]) washes in a customized microscope stage), the CODEX tech- and simultaneous placement of those cells in the tissue architec- nique could enable deep studies of various tissue models even ture. Cell-interaction analysis with CODEX recapitulated known with limited resources and instrumentation. features of splenic tissue architecture and revealed that most The unique set of algorithms described here successfully splenic cell types were frequently in homotypic interactions— identified individual cells in the crowded environment of lymphoid which might underscore a novel driving principle of lymphoid tis- tissue by relying on both the information from nuclear and the sue architecture. Further, an important principle was derived from membrane staining. An accurate quantification of single-cell the data-driven i-niche analysis; i.e., we have established that expression data was obtained directly from the images by certain markers, such as B220, CD79b, or CD27, exhibited signif- creating a special algorithm for positional spill compensation. As icant changes in expression levels depending on the tissue of today, this algorithm is only applicable to uniformly distributed context in which the cells reside. As clearly observed in experi- surface markers. Future changes might be required to accommo- ments that drove Figures 4 and 6, cell populations that would date markers that follow a different distribution, i.e., localized to otherwise be thought of as ‘‘broadly’’ expressing a given marker lipid rafts or immune synapses. Nevertheless, the use of this algo- set(Figure 4H), in fact, were composed of multiple subphenotypes rithm enabled us to extract FACS-like data from tissue imaging that correlated with the i-niche identity. In other words, what im- and leveraged the automated phenotype mapping framework munologists previously thought of as a single cell type could be previously developed for CyTOF and multicolor FACS. subdivided into more subtle cell subsets that are defined by the Performance of CODEX on tissue sections was validated in neighborhood in which they reside. We leave open the question analysis of spleen sections of normal and lupus afflicted mice of whether the cells with different properties are attracted to a Cell 174, 968–981, August 9, 2018 979 d EXPERIMENTAL METHOD DETAILS set of neighbor cells, or a given expression level of markers at- B Oligonucleotide sequences tracts the neighbors, or some dialectics thereof. What is clear, B Primer dependent panels to extend the multiplexing however, is that there are more subtle phenotypes in tissues capacity of CODEX than previously assumed, and that future developments of the B CyTOF CODEX comparison technology and the algorithms will shed more light on these B Antibody conjugation, staining and CODEX rendering phenomena. B Primer dependent panels CODEX enabled a quantitative description of autoimmune- B Imaging related changes in the splenic tissue architecture. Among hall- d QUANTIFICATION AND STATISTICAL ANALYSIS marks of MRL/lpr progression were dissipation of marginal B Initial image analysis and segmentation zone, disintegration of PALS, invasion of red pulp with erythro- B Spatial spillover compensation blasts, and the infiltration of mixed-identity B220 DN T cells, B Cell type definition which, interestingly, localize in a niche in between PALS and B Cell interaction analysis the B cell zone and in the marginal zone. A contact-dependent B i-niche analysis effect of B220 DN T cell on CD4 T cells reflected in increased B Neural network training and data analysis levels of activation marker CD27 was observed. An account of d DATA AND SOFTWARE AVAILABILITY statistically significant differences in frequency and strength of pairwise cell-type contacts was created. From these observa- SUPPLEMENTAL INFORMATION tions and their quantitative analysis, we concluded that it is largely the change in cell numbers rather than in cellular interac- SupplementalInformation includesseven figures, fourtables,and two videosand tion strength (estimated from ratio of observed to expected can be found with this article online at https://doi.org/10.1016/j.cell.2018.07.010. probability of interaction) that is involved in reorganization of spleen during transition from norm to autoimmunity. We show ACKNOWLEDGMENTS how i-niche statistics can be used to account for the list of dis- ease-driven changes in sub-splenic anatomy. We also show We thank A. Trejo and A. Jager for technical assistance. We are grateful to Peter Jackson for advice and access to the microscope used for this work. that disease-associated areas of the tissue can be identified This work was further supported by the US NIH, grants and sub awards: independently of the image segmentation, by applying a convo- U19AI057229, 1U19AI100627, R01CA184968, 1R33CA183654-01, lutional neural network to the multidimensional image data, even R33CA183692, 1R01GM10983601, 1R21CA183660, 1R01NS08953301, after training on just one sample. 5UH2AR067676, 1R01CA19665701, R01HL120724, and R01HL128173, Bill Recent advances in genomics suggest that, despite the vast- and Melinda Gates foundation grant OPP1113682, Department of Defense ness of a genetic repertoire, there exist only a limited number of (CDMRP) grant W81XWH-14-1-0180, Northrop-Grumman Corporation agree- cellular states with a concomitantly limited gene expression ment 7500108142, and NIH grant to Stanford Center of Cancer Systems Biology U54-CA209971. G.P.N. is supported by the Rachford & Carlotta A. pattern. These countable, limited patterns are reflected in the Harris Endowed Chair. expression of surface-marker phenotypes recognizable as cell types. It is therefore reasonable to suggest that cell-to-cell AUTHOR CONTRIBUTIONS interactions should be limited as well and falling into repeated patterns. By this token, the data collected in this study lays Conceptualization, Y.G., N.S., and G.P.N.; Methodology, Y.G., N.S., and the foundation for a pan-cellular reference database defining J.K.-D.; Software, N.S. and S. Bhate; Validation, Y.G., N.S., J.K.-D., and G.V.; Formal Analysis, Y.G. and N.S.; Investigation, Y.G. and N.S.; Resources, cellular types not only by identities of proteins expressed, but Y.G., N.S., G.V., S. Black, and M.H.; Data Curation, Y.G. and N.S.; Writing – also by definitions for specific cell-to-cell interactions. We Original Draft, Y.G. and N.S.; Writing – Review and Editing, Y.G., N.S., and performed deep characterization here for normal and diseased G.P.N.; Visualization, Y.G. and N.S.; Supervision, Y.G. and G.P.N.; Project tissue from such a perspective of cell-cell arrangements and Administration, Y.G. and G.P.N.; Funding Acquisition, G.P.N. present here, for the research community, a large (700,000 cells) public dataset encompassing segmentation, quantifica- DECLARATION OF INTERESTS tion, and, most uniquely, spatial data from normal and disease- afflicted spleens (http://welikesharingdata.blob.core.windows. G.P.N., Y.G., and N.S. all have ownership in or are consultants of Akoya Bio, Inc., and members of its scientific advisory board. J.K.-D. was at Stanford net/forshare/index.html). Further analysis of data could enable during the study and currently is an employee of Akoya Bio. Akoya makes re- advances in understanding of clinically relevant cell interactions agents/instruments that are dependent upon licenses from Stanford Univer- in immune tissues as well as development of computational al- sity. Stanford University has been granted a US patent 9909167 that covers gorithms for tissue cytometry and digital pathology. some aspects of the technology described in the paper. Received: September 22, 2017 STAR+METHODS Revised: February 5, 2018 Accepted: July 3, 2018 Detailed methods are provided in the online version of this paper Published: August 2, 2018 and include the following: REFERENCES d KEY RESOURCES TABLE d CONTACT FOR REAGENT AND RESOURCE SHARING Bendall, S.C., Simonds, E.F., Qiu, P., Amir, A.D., Krutzik, P.O., Finck, R., d EXPERIMENTAL MODEL AND SUBJECT DETAILS Bruggner, R.V., Melamed, R., Trejo, A., Ornatsky, O.I., et al. (2011). Single-cell 980 Cell 174, 968–981, August 9, 2018 mass cytometry of differential immune and drug responses across a human he- Parslow, A., Cardona, A., and Bryson-Richardson, R.J. (2014). Sample drift matopoietic continuum. Science 332, 687–696. correction following 4D confocal time-lapse imaging. J. Vis. Exp. Published online April 12, 2014. https://doi.org/10.3791/51086. Burrell, B.E., Warren, K.J., Nakayama, Y., Iwami, D., Brinkman, C.C., and Preibisch, S., Saalfeld, S., and Tomancak, P. (2009). Globally optimal stitching Bromberg, J.S. (2015). Lymph node stromal fiber ER-TR7 modulates CD4+ of tiled 3D microscopic image acquisitions. Bioinformatics 25, 1463–1465. T cell lymph node trafficking and transplant tolerance. Transplantation 99, Robinette, M.L., Fuchs, A., Cortez, V.S., Lee, J.S., Wang, Y., Durum, S.K., Gil- 1119–1125. fillan, S., and Colonna, M.; Immunological Genome Consortium (2015). Tran- Chattopadhyay, P.K., and Roederer, M. (2012). Cytometry: Today’s technol- scriptional programs define molecular characteristics of innate lymphoid cell ogy and tomorrow’s horizons. Methods 57, 251–258. classes and subsets. Nat. Immunol. 16, 306–317. Gabriel, K.R., and Sokal, R.R. (1969). A new statistical approach to geographic Rubtsov, A.V., Rubtsova, K., Fischer, A., Meehan, R.T., Gillis, J.Z., Kappler, variation analysis. Syst. Zool. 18, 259. J.W., and Marrack, P. (2011). Toll-like receptor 7 (TLR7)-driven accumulation of a novel CD11c B-cell population is important for the development of auto- Gerdes, M.J., Sevinsky, C.J., Sood, A., Adak, S., Bello, M.O., Bordwell, A., immunity. Blood 118, 1305–1315. Can, A., Corwin, A., Dinn, S., Filkins, R.J., et al. (2013). Highly multiplexed sin- gle-cell analysis of formalin-fixed, paraffin-embedded cancer tissue. Proc. Rubtsova, K., Rubtsov, A.V., Thurman, J.M., Mennona, J.M., Kappler, J.W., Natl. Acad. Sci. USA 110, 11982–11987. and Marrack, P. (2017). B cells expressing the transcription factor T-bet drive lupus-like autoimmunity. J. Clin. Invest. 127, 1392–1404. Gerner, M.Y., Kastenmuller, W., Ifrim, I., Kabat, J., and Germain, R.N. (2012). Samusik, N., Good, Z., Spitzer, M.H., Davis, K.L., and Nolan, G.P. (2016). Auto- Histo-cytometry: A method for highly multiplex quantitative tissue imaging mated mapping of phenotype space with single-cell data. Nat. Methods 13, analysis applied to dendritic cell subset microanatomy in lymph nodes. Immu- 493–496. nity 37, 364–376. Schubert, W., Bonnekoh, B., Pommer, A.J., Philipsen, L., Bo¨ ckelmann, R., Ma- Jacobson, B.A., Panka, D.J., Nguyen, K.A., Erikson, J., Abbas, A.K., and lykh, Y., Gollnick, H., Friedenberger, M., Bode, M., and Dress, A.W. (2006). Marshak-Rothstein, A. (1995). Anatomy of autoantibody production: Dominant Analyzing proteome topology and function by automated multidimensional localization of antibody-producing cells to T cell zones in Fas-deficient mice. fluorescence microscopy. Nat. Biotechnol. 24, 1270–1278. Immunity 3, 509–519. Socolovsky, M., Murrell, M., Liu, Y., Pop, R., Porpiglia, E., and Levchenko, A. Kanauchi, H., Furukawa, F., and Imamura, S. (1991). Characterization of cuta- (2007). Negative autoregulation by FAS mediates robust fetal erythropoiesis. neous infiltrates in MRL/lpr mice monitored from onset to the full development PLoS Biol. 5, e252. of lupus erythematosus-like skin lesions. J. Invest. Dermatol. 96, 478–483. Spitzer, M.H., Gherardini, P.F., Fragiadakis, G.K., Bhattacharya, N., Yuan, Koh, D.R., Ho, A., Rahemtulla, A., Fung-Leung, W.P., Griesser, H., and Mak, R.T., Hotson, A.N., Finck, R., Carmi, Y., Zunder, E.R., Fantl, W.J., et al. T.W. (1995). Murine lupus in MRL/lpr mice lacking CD4 or CD8 T cells. Eur. (2015). IMMUNOLOGY. An interactive reference framework for modeling a dy- J. Immunol. 25, 2558–2562. namic immune system. Science 349, 1259425. Lieberum, B., and Hartmann, K.U. (1988). Successive changes of the cellular Spitzer, M.H., Carmi, Y., Reticker-Flynn, N.E., Kwek, S.S., Madhireddy, D., composition in lymphoid organs of MRL-Mp/lpr-lpr mice during the develop- Martins, M.M., Gherardini, P.F., Prestwood, T.R., Chabon, J., Bendall, S.C., ment of lymphoproliferative disease as investigated in cryosections. Clin. Im- et al. (2017). Systemic immunity is required for effective cancer immuno- munol. Immunopathol. 46, 421–431. therapy. Cell 168, 487–502. Cell 174, 968–981, August 9, 2018 981 STAR+METHODS KEY RESOURCES TABLE REAGENT or RESOURCE SOURCE IDENTIFIER Antibodies B220 (Ra3-6B2) BD Cat#: 553084; RRID: AB_394614 CD106 (429(MVCAM.A)) BD Cat#: 553330; RRID: AB_394786 CD11b (M1/70) BD Cat#: 553308; RRID: AB_394772 CD11c (HL3) BD Cat#: 117302; RRID: AB_313771 CD16/32 (2.4G2) BD Cat#: 553142; RRID: AB_394657 CD169 (MCA947G) BioRad Cat#: MCA947G; RRID: AB_322322 CD19 (1D3) BD Cat#: 553783; RRID: AB_395047 CD21/35 (7G6) BD Cat#: 553817; RRID: AB_395069 CD27 (LG.3A10) BD Cat#: 124202; RRID: AB_123645 CD3 (17A2) BD Cat#: 555273; RRID: AB_395697 CD31 (MEC13.3) BD Cat#: 553370; RRID: AB_394816 CD35 (8C12) BD Cat#: 558768; RRID: AB_397114 CD4 (RM4-5) BD Cat#: 553043; RRID: AB_394579 CD44 (IM7) BD Cat#: 553131; RRID: AB_394646 CD45 (30-F11) BD Cat#: 553076; RRID: AB_394606 CD5 (53-7.3) BioLegend Cat#: 100602; RRID: AB_312731 CD71 (C2F2) BD Cat#: 553264; RRID: AB_394742 CD79b (HM79b) BioLegend Cat#: 132802; RRID: AB_207588 CD8a (53-6.7) BD Cat#: 553027; RRID: AB_394565 CD90 (G7) BioLegend Cat#: 105202; RRID: AB_313169 ERTR7 (ERTR7) BioRad Cat#: MCA2402; RRID: AB_915429 F4/80 (T45-2342) BD Cat#: 565409 IgD (11-26c.2a) BD Cat#: 553438; RRID: AB_394858 IgM (II/41) BD Cat#: 553405; RRID: AB_394842 Ly6C (HK1.4) Biolegend Cat#: 128001; RRID: AB_113421 Ly6G (1A8) BD Cat#: 551459; RRID: AB_394206 MHCII (M5/114.15.2) BD Cat#: 556999; RRID: AB_396545 Ter119 (TER-119) Biolegend Cat#: 116202; RRID: AB_313703 For composition of oligo conjugated and metal-labeled This paper N/A abs see Table S1 and antibody conjugation part of STAR Methods Biological Samples Spleens from MRL/lpr mice This paper N/A Spleens from BALBc mice This paper N/A Isolated mouse splenocytes This paper N/A Chemicals, Peptides, and Recombinant Proteins dUTP-ss-Cy5; dCTP-ss-Cy3 Jena Custom order Exo- fragment of Klenow DNA polymerase NEB Cat#: M0212L TCEP SIGMA-Aldrich Cat#: C4706-10G Iodoacetamide SIGMA-Aldrich Cat#: I1149-5G Experimental Models: Organisms/Strains MRL/lpr mice Jackson Laboratory IMSR Cat# JAX:000485, RRID:IMSR_JAX:000485 BALBc Jackson Laboratory IMSR Cat# JAX:000651, RRID:IMSR_JAX:000651 (Continued on next page) e1 Cell 174, 968–981.e1–e6, August 9, 2018 Continued REAGENT or RESOURCE SOURCE IDENTIFIER Oligonucleotides (see Table S3) N/A N/A Software and Algorithms CODEX Toolkit This paper https://github.com/nolanlab/CODEX VORTEX environment for X-shift clustering This paper https://github.com/nolanlab/VORTEX Microvolution software for image deconvolution Microvolution http://www.microvolution.com Deposited Data Dataset with data tables can be found at Mendeley https://data.mendeley.com/datasets/zjnpwh8m5b/1 Other Primary image data and additional info can be found at Online repository https://welikesharingdata.blob.core.windows.net/ forshare/index.html CONTACT FOR REAGENT AND RESOURCE SHARING Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Prof. Dr. Garry Nolan (gnolan@stanford.edu). EXPERIMENTAL MODEL AND SUBJECT DETAILS 9 months old female MRL/lpr (IMSR Cat# JAX:000485, RRID:IMSR_JAX:000485) (chosen to represent lupus disease at a pronounced splenomegaly stage) and age/sex matched control BALBc mice (IMSR Cat# JAX:000651, RRID:IMSR_JAX:000651) purchased from Jackson Laboratory were used for the study. All animal studies were done in compliance with ethical regulations and procedures set in the Stanford Administrative Panel on Laboratory Animal Care Protocol 15986. In coherence with the primarily technical purpose of the study no animal cohort randomization or investigator blinding to group allocation was performed. EXPERIMENTAL METHOD DETAILS Oligonucleotide sequences Single base extension during CODEX can be achieved by either a ‘‘missing base’’ approach (Figure 1A) or a ‘‘reversible terminator’’ method (see Video S1, part 2). In the case of the ‘‘missing base’’ approach, which was chosen for the experiments outlined in this paper, the top strand of the double-stranded oligonucleotide is covalently bound to the capture agent (in this case, an antibody) and the bottom strand is annealed through hybridization to the top strand. All antibodies contain the same top strand (5 -ATAG CAGTCCAGCCGAACGGTAGCATCTTGCAGAA-3 ) and different bottom strands. The sequence of the bottom strands contains a common region that hybridizes to the top strand (.TTCTGCAAGATGC 0 0 TACCGTTCGGCTGGAddC-3 ) as well as a 5 variable sequence region that serves as the indexing region. As shown in Figure 1A, the overhanging 5 end of the lower strand of the double-stranded oligonucleotide tag (which forms the overhang) is of the general 0 0 0 0 formula 5 -[C/T] [A/G][5 -C /T -3 ] TTCTGCAAGATGCTACCGTTCGGCTGGAddC-3 The first block a short 5-nt stretch of 5 1-4 1-4 n- random C/T composition designed to increase the polymerase residence on the DNA duplex. The second block is a single nucleotide (either G or A) that allows for incorporation of a labeled dNTPs (dU-ss-Cy5 or dC-ss-Cy3, respectively). The third block is the ‘‘index- ing barcode‘‘ that consists of n random-length homopolymer stretches (1-4 nucleobases each) of alternating ‘‘indexing’’ nucleobases dC and dT that serve as a template for extension of the top oligo with unlabeled nucleotides (dATP and dGTP). Here, n specifies the number or extension cycles after which the fluorescent nucleobase will be incorporated into the duplex. Examples of CODEX index- ing barcodes are CCCTCC for n = 3 and CCCTCCTTTCTT for n = 6. The purpose of having the homopolymer stretches of variable length (e.g., CCCTCCTTTCTT) rather than single base (e.g., CTCTCT) is to increase the polymerase extension specificity and prevent misalignment of upper and lower strands of double-stranded oligonucleotide tags. All oligonucleotide sequences used in this study can be found in Table S1. Primer dependent panels to extend the multiplexing capacity of CODEX CODEX operates using an indexed polymerization step that enables precise incorporation of fluorophores into oligonucleotide-Ab conjugates at predetermined cycles. Although consistent performance of a model antigen (CD45) was observed across 15 cycles of CODEX (Figures S1A–S1F), a gradual accumulation of polymerization errors during each cycle could potentially result in Cell 174, 968–981.e1–e6, August 9, 2018 e2 non-cognate rendering, and thus diminished and/or non-specific signals at later index cycles. In addition, the use of long single- stranded oligonucleotides that would enable indexing beyond 15 rounds might be problematic due to non-specific binding events to tissues under study. For the polymerization event to initiate, a 3 hydroxyl is required. Thus, we reasoned that dedicated primers (each containing a distinct initiating sequence with a 3 hydroxyl) could be used to activate distinct subpanels of antibodies (Figure S2A). This would allow design of antibody panels exceeding 30 markers into subpanels, each with a subpanel-specific activation sequence designed 5 to the indexing region. In this design, the antibody attachment linker is terminated with ddC, such that the extension is only possible after a hybridization of a hydroxyl-containing panel-specific activation primer. The feasibility of such multipanel CODEX design and the robustness of CODEX protocol after many cycles and its independence of staining from the cycle number were tested in a model experiment. A 22-color panel of antibodies (11 cycles) conjugated to a termi- st nd rd nated top oligonucleotide, was hybridized with lower oligonucleotides of 1 ,2 , and 3 panels (Figure S2B). Thus, every antigen is detected thrice by the same antibody conjugated to oligonucleotides of 3 different panels. Each panel can only be rendered after annealing of a panel-specific activator oligonucleotide. The staining was rendered in 36 cycles (11 detection cycles + 1 blank no-anti- body cycle per activator oligo) of CODEX with additional activator oligonucleotide hybridization step between each of the 3 panels. st th th The signal for same antibody detected at different cycles (e.g., 1 ,13 , and 24 ) was consistent across the three panels (Figure S2C). This panel-activator design extends CODEX to a theoretically unlimited multiplexing capacity, bounded only by the speed and res- olution of the imaging process itself and the time required for each imaging cycle. CyTOF CODEX comparison Cell preparation and staining by metal tagged antibody for CyTOF analysis was performed as described before (Spitzer et al., 2017). TM Mass cytometry was performed on a CyTOF 2 mass cytometer (Fluidigm) equilibrated with ddH2O. For CODEX analysis, isolated spleen cells were stained a panel of antibodies conjugated to indexing oligonucleotides. Samples were fixed to a coverslip (Figure 2A) and imaged over 12 cycles of CODEX protocol. Images were segmented using the in situ cytometry software toolkit developed for this study (see Figure 2A for exemplary segmentation of the cell spread), and the staining of individual cells across the indexing cycles was quantified. Segmentation data was converted into flow cytometry standard (FCS) format and analyzed using the conventional flow cytometry analysis software Cytobank. Antibody conjugation, staining and CODEX rendering Detailed stepwise CODEX protocols can be found online (see Key Resources Table and Data and Software Availability section below). For full list of antibody clones and vendors see Table S1. Custom manufactured microfluidic setup (Figures S7A–S7C) was used to automate CODEX solution exchange and image acquisition. Instrument and blueprints and control software are available upon request. Primer dependent panels Rendering of antibodies with spacers followed the same procedure as the standard CODEX protocol with the exception of the following differences. Before proceeding to rendering next spacer dependent panel, the stained cells were incubated with a spacer oligonucleotide (1 mM final concentration in buffer 405) at room temperature for 10 minutes. Cells were washed 4X with buffer 4 and rendering proceeded as usual. To initiate each additional spacer set, the spacer incubation step was repeated using corresponding spacer samples. Imaging Images were collected using a Keyence BZ-X710 fluorescent microscope configured with 3 fluorescent channels (FITC, Cy3, Cy5) and equipped with Nikon PlanFluor 40x NA 1.3 oil immersion lens. Imaging and washes were iteratively performed automatically us- ing a specially developed fluidics setup (Figures S7A–S7C). Images were subject to deconvolution using Microvolution software (http://www.microvolution.com/). The staining patterns of 28 DNA-conjugated antibodies were acquired over 14 cycles of CODEX imaging and overlaid with 2 additional fluorescent antibodies, CD45-FITC and NKp46-Pacific Blue and a DRAQ5 nuclear stain (Figure 3A and low-resolution views in Video S2). Each tissue was imaged with a 40x oil immersion objective in a 7x9 tiled acquisition at 1386x1008 pixels per tile and 188 nm/pixel resolution and 11 z-planes per tile (axial resolution 900 nm). Images were subjected to deconvolution to remove out-of-focus light. After drift-compensation and stitching, we obtained a total of 9 images (one per tissue) with x = 9702 y = 9072 z = 11 dimensions, each consisting of 31 channels (30 antibodies and 1 nuclear stain). QUANTIFICATION AND STATISTICAL ANALYSIS Initial image analysis and segmentation For each imaging field analyzed by CODEX multidimensional staining multi-color z stacks collected during individual cycles were aligned against reference channel (CD45) by 3D drift compensation (Parslow et al., 2014). If necessary individual fields covering large tiled areas were ‘‘stitched’’ using dedicated ImageJ plugin (Preibisch et al., 2009). For the 22-channel experiment on dissociated cells attached to coverslip (Figure 2) images corresponding to the best focal plane of vertical image stacks collected at each acquisition e3 Cell 174, 968–981.e1–e6, August 9, 2018 step of CODEX were chosen for quantification. For the 31-channel main experiments on mouse spleen sections, the segmentation was performed on the whole image stack using a volumetric segmentation algorithm described below. For this study, we purposefully developed a 3D image segmentation that combines information from the nuclear staining and a ubiquitous membrane marker (in this case CD45) to define single-cell boundaries in crowded images such as lymphoid tissues. This algorithm inverts the membrane image and multiplies it with the nuclear image, creating a synthetic image with enhanced contrast between neighboring nuclei. This image is subject to low-pass FFT filtering and an then individual cell objects (collections of voxels) are identified using a gradient-tracing water- shed algorithm. Per-cell intensities were quantified by integrating the intensity of each channel within a given cell object and divided by the region size in voxels. We benchmarked the segmentation algorithm against a dataset of BALBc mouse spleen images with expert hand-labeled nuclei in and we found the algorithm was able to correctly identify 87.25% ± 2.89% cells, of which 89.88% ± 1.12% were singlets (one-to-one correspondence between a hand-labeled cell and a segmented object) (Figures S1H–S1J). For each segmented object (i.e., cell) a marker expression profile, as well as the identities of the nearby neighbors were recorded using Delaunay triangulation (https://data. mendeley.com/datasets/zjnpwh8m5b/1). Spatial spillover compensation Accurate segmentation per se is not sufficient to obtain high-quality estimate of single-cell expression from images. The reason for that is that in lymphoid tissue the cells are so tightly adjacent that their membrane signals can partially overlap, resulting in blending of signals between neighboring cells, the phenomenon termed spatial spillover. In order to compensate for that, we estimated the cell- to-cell signal spill coefficients based on the fraction of shared boundary between each pair of cell objects, resulting in a banded matrix (most cells don’t have any shared boundaries). To compensate the cell-to-cell spill, the raw intensity vector was multiplied by the inverse spill matrix (Figure 2D). Besides the spatial signal spillover, there are other factors that add to artifactual cell-like objects: debris misidentified as cells, dou- blets (two adjacent cells merged together) as well as autofluorescent objects, both of which can lead to spurious as double-positive signals on the biaxial scatterplots. By analogy to how debris and doublets are eliminated from FACS data by applying special ‘Singlet’ gates to SSC-FSC parameters, we devised a ‘cleanup’ gating strategy based on several quality control parameters: nuclear stain density (nuclear signal divided by cell size), profile homogeneity (relative variance of signal from cycle to cycle), background staining on blank cycles and, finally, nuclear signal and cell size (Figure S1K). We found that applying those filtering gates had a synergetic effect with the compensation, reducing the frequency of spurious double-positive cell signals by approximately an additional factor of 2 (Figure S1L, e.g., compare fraction of CD5/IgD double positive cells in Ungated-Compensated and Post Cleanup gated – Compensated in (L)). Cell type definition The 9-spleen dataset was subject segmentation, quantification, compensation and cleanup gating, as described above, yielding a total of 734101 30-dimensional single-cell protein marker expression profiles (Figure 3C, https://data.mendeley.com/datasets/ zjnpwh8m5b/1). The segmented CODEX data was subject to automated phenotype mapping algorithm X-shift that was previously developed and validated on CyTOF data (Samusik et al., 2016)(Figure 3C). 58 phenotypic clusters inferred by X-shift clustering were manually annotated (Figures 3C and 3D and Table S2) based on the 30-color marker expression profile and thorough visual inspec- tion of the representative image samples (Figures S5A–S5J, more examples of ‘‘cell passports’’ can be found in associated online repository - see STAR Methods). Some clusters were found to originate from imaging artifacts such as dust and tissue sectioning defects. That reduced the overall number of cell-like objects to 707466. Each cluster was assigned to one of 27 broadly defined sin- gle-cell phenotypic groups (cell types), which in some cases could be clearly matched to major immune cell types and in others were named according to expression of distinguishing surface markers (see cluster annotation and cell counts in Table S2). Cell interaction analysis To define for each cell the neighbors of the first (immediate) tier of proximity Delaunay graph was computed for the dataset (https:// data.mendeley.com/datasets/zjnpwh8m5b/1). The odds ratio of co-occurrence of cell type A and cell type B was estimated as the observed frequency of co-occurrence (mean of the beta-distribution, with parameter alpha = number of edges connecting cell types A and B and parameter beta = total number of edges minus number of edges connecting A-B) divided by the theoretical frequency of co-occurrence (total frequency of edges incident to type A multiplied by the total frequency of edges incident to type B) see Table S3. The odds ratios are represented in heatmaps on Figure 3G, with a range of values from less than 1 to more than 1 mean- ing that two cell types are, respectively, less or more likely to co-occur than expected by chance. The significance of the difference from zero was tested using binomial distribution (probability of getting an observed number of interactions between A and B (successes) among the total number of registered interactions (number of trials) given the theoretical probability of A-B interaction (probability of success)). The significance of change of interaction frequencies or log-odds ratios were computed between BALB/c and Stage 1 (early) MRL using pairwise t test. However, the same procedure could not be applied to testing BALB/c versus MRL/lpr Stages 2 and 3 because of high sample-specific variation in those more advanced disease stages. Therefore we scored computed the deviation of those Stage2/3 values from BALB/c using c statistics because it does not require Stage 2/3 samples to have a common mean. Cell 174, 968–981.e1–e6, August 9, 2018 e4 The P values were subject to FDR correction using Benjamini–Hochberg procedure. Interactions that were considered significant for FDR q-value < 0.05 or > 0.95 (Table S3). In order to estimate the overall deviation of either interaction frequency matrices or log-odds ratio matrices, the matrices were sub- ject to z-transformation based on the mean and the SDs of the BALB/c samples, and then c statistics was computed as square root of the sum of squares of all elements of the z-score transformed matrices (Figure 5F). i-niche analysis The i-niche analysis was performed based on 2-dimensional Delaunay triangulation of the cell center coordinates. Delaunay triangu- lation and the related concepts of Voronoi Tesselation and Gabriel graphs were previously applied in eco-geographical analysis of species distribution (Gabriel and Sokal, 1969) and therefore were deemed as equally applicable to the analysis of tissue organization on the single-cell level. The i-niche is defined as a set of first-order Delaunay neighbors of the given ‘index cell’, i.e., the i-niche cells are the ones that are directly connected to the i-cell with edges in the Delaunay triangulation of cell centers. We distinguish i-niche from the more formal understanding of ‘‘niche,’’ which is often used in stem cell literature and where numbers of cells in the niche and their placement within the niche is undefined. In our definition, we allow the central cell to be of any type and are counting the cell types present in the ring. This flexible definition allows for multi-cellular interactions around a central cell to define the biology of that cell (and vice versa). Computationally, the i-niche window slides from cell to cell, considering each set of adjoining cells— and therefore allows consideration of the constituencies of different central cell types that might populate a given i-niche. We under- stand that our current definition is arbitrary and could be extended to include other specific cell arrangements—including, though beyond the scope of the current work, a 3D sphere of cells contacting the index cell. Neural network training and data analysis Preprocessing Image stacks were maximum-intensity projected following deconvolution. Data was quantile normalized to 4 levels (0, 0.25, 0.5 and 0.75 quantiles). A baseline model was able to distinguish models without this discretization and normalization, suggesting strain-spe- cific differences in antibody staining intensity. Training and cross validation split Four spleen samples (two BALB/c and two MRL/lpr) were chosen as training samples. The remaining five spleens tissue samples (one BALB/c and four MRL/lpr) were used for testing the trained model. For cross-validation, different combinations of spleens were allo- cated to training and test sets. During training, 224x224 images were randomly extracted from the training tissue samples, at 1x, 0.5x and 2x zoom. At 1x zoom, there would be 6804 non-overlapping image patches in the training dataset. The trained models were tested on 4500 patches, at 1x zoom. Hyperparameters were manually tuned on 500 randomly selected images from the testing spleens. The Adam optimizer was used for training with an initial learning rate of 0.0001. Baseline model. A logistic regression model was trained by averaging marker intensities across the image. L2 regularization was used for weights. Neural network architecture. To avoid the learning of trivial sample-specific staining variation, data were quantile normalized sam- ple-wise and each marker was discretized to four levels. Since disease-specific hallmarks could potentially be present at multiple scales, the training data for our neural network was extracted at multiple levels of magnification. A simple regularized logistic regres- sion model that considered only average marker expression and did not incorporate spatial information was unable to successfully distinguish patches normal and MRL/lpr spleens, whereas the trained neural network model consistently achieved a 90% precision of classification of image patches during cross-validation. A fully convolutional network architecture was used, with the following layers. To generate a prediction for an entire image patch, a global max-pooling layer was used. 1. Conv3 60 2. Conv3 120 3. Conv3 64 4. Batch Norm 5. Conv3, 64 6. Max pooling 2x2 7. Conv3, 128 8. Conv3, 128 9. Max pooling 2x2 10. Conv3, 256, 11. Conv3,256 12. Conv3,256 13. Max pooling 2x2 14. Conv3,512 15. Conv3,512 e5 Cell 174, 968–981.e1–e6, August 9, 2018 16. Conv3,512 17. Conv1,256 18. Conv1,64 19. Conv1,1 20. Global max pooling 21. Sigmoid Weights for layers 5-16 were initialized from the VGG-16 pretrained model. The model was trained with cross-entropy loss. Regularization. L2 regularization (0.1) was used for network weights. L1 regularization was applied to the feature map output after layer 19 to encourage sparse activations Whole sample activations for test set. Since the network was fully convolutional, it could be applied to images of any dimension. The network was applied to entire fields of view individually. The activation maps were obtained as the output after layer 21. Aligning cell type information. Each cell was assigned the MRL/lpr score of the corresponding pixel in the image. Enrichment and neighborhood analysis. FDR controlled chi-square tests of proportions were carried out to determine enrichment of specific cell types in the top 10% of cells by MRL/lpr score. For neighborhood analysis of dendritic cells, the composition of the neighborhoods (cell centers within 30 pixels) of the top 300 cells (by MRL/lpr score) were compared to the composition of the neighborhoods of the bottom 300 cells. Only cells with positive neural network assigned MRL/lpr score, in MRL/lpr regions, were considered for this analysis. DATA AND SOFTWARE AVAILABILITY Software used in the paper for parsing image data can be obtained at: https://github.com/nolanlab/CODEX Data tables can be downloaded from Mendeley: https://data.mendeley.com/datasets/zjnpwh8m5b/1 All primary image data, high resolution focused montages, complete single cell data tables and various additional information can be obtained at: http://welikesharingdata.blob.core.windows.net/forshare/index.html Flow formatted segmented data can obtained from online repository page (link above) or from Cytobank: CODEX on spreads of isolated mouse splenocytes (Figure 2B): https://community.cytobank.org/cytobank/experiments/69534 https://flowrepository.org/experiments/1686 CyTOF on isolated splenocytes (Figure 2B): https://community.cytobank.org/cytobank/experiments/69533 https://flowrepository.org/experiments/1687 CODEX on BALBc spleen tissue sections (Figure 2E): https://community.cytobank.org/cytobank/experiments/69889 https://flowrepository.org/experiments/1688 Cell 174, 968–981.e1–e6, August 9, 2018 e6 Supplemental Figures Figure S1. Benchmarking CODEX, Related to Figures 1 and 2 (A) Experimental scheme for mimicking the tissue with 30 distinct cell types. (B) Montage of a fragment of imaging field of the 15 cycles of CODEX used to render the mix of 30 barcoded spleens – first cycle top left last cycle bottom right. (C) Heatmap (cycles in columns, cells in rows) showing mean fluorescence per cell membrane for each cell per in each of the 15 CODEX cycles performed on cells of 30 barcoded spleens. Odd columns correspond to imaging after labeled base incorporation. Even columns correspond to imaging after inactivation of staining by TCEP. (D) Time-lapse profile of median intensity per cell membrane for individual cells marked by white arrows on (B). (E) Average intensity of CD45 antigen expression in ‘‘positive’’ (blue columns) and ‘‘negative’’ (red columns) cells in 15 CODEX cycles of the experiment. (Similar results were obtained for Cy3-positive populations – data not shown). Linear regression was performed to indicate trends in accumulation of background and signal decline associated with cycle number. (F) Table summarizing CODEX performance stats. Average signal to noise ratio was estimated from ratio of average signal of all positive cells across all cycles to the signal of all ‘‘negative’’ cells across all cycles. Efficiency of fluorophore removal was estimated from average ratio of ([signal after TCEP in cycle N]- [signal after TCEP in cycle (N-1)])/[signal in cycle N] for cells positive in cycle N across all cycles. Average expected signal deterioration was estimated using the trendline equation from (E). Average background accumulation was estimated by fitting linear trendline into the per cycle ratio of average background to average signal (not shown). (G) Image quantification approach used on CODEX data from (A): best focal planes of CODEX stacks were segmented by Cell Profiler. To account for local background the value corresponding to difference between the mean intensity value inside ‘‘cell membrane’’ object (left panel) and the mean intensity inside the (legend continued on next page) external ring object (right panel) was chosen as a representation of the intensity of the antibody signal. In all other experiments custom (see STAR Methods) segmantation developed in this study was used. (H) Sample 500x500 px regions from two samples (BALBc-3 and MRL-4) showing hand-labeled cell centers (yellow crosshairs) and cell outlines detected by the segmentation algorithm (randomly colored). (I) Comparison between the hand-labeled cell identification and algorithm-based algorithm identification, expressed in 3 measures: %Nuclei found (how many of the hand-labeled nuclei centers ended up inside the segmented regions), % Singlets (how many of the cell regions with at least one hand-labeled nuclei center contained exactly one cell center) and %Unlabelled regions (how many segmented regions did not contain a hand-labeled cell center). (J) Summary statistics comparing the segmentation quality between BALBc and MRL/lpr samples. (K) Three step cleanup gating strategy based on 1) stain density (nuclear signal divided by cell size) and profile homogeneity (relative variance of signal from cycle to cycle), 2) removing objects with high background by gating on the signal accrued in ‘‘blank’’(no stain) cycles 3) constraining the cell size. (L) Percentage of artifactual double-positive cells in CODEX data from sample BALBc-2 (as seen in the upper right quadrant biaxial flow style plots of mutually exclusive lineage markers IgD and CD5) depending on gating and spill compensation. Figure S2. Expanding the Multiplexing Limit of CODEX by ‘‘Panels and Activators’’ Design, Related to Figure 1A and STAR Methods (A) Diagram of ‘‘multipanel’’/’’activator oligo’’ CODEX approach. The list of antibodies can be divided in sets such that number of antibodies in each individual set does not exceed the capacity of the multiplexing protocol to render staining without significant signal loss (e.g.30). Each such set of antibodies will be conjugated to ‘‘terminated’’ (the last 3 base is dideoxy- or propyl- modified) upper strand oligonucleotide of the same sequence as in the original version of the ‘‘missing base’’ approach. The lower strand oligonucleotides will incorporate an additional set-specific region, which will serve as a landing spot for the dedicated primer oligo which is to be on-slide hybridized to the particular subset of the total plurality of the antibodies at the time when they are to be rendered. This approach prevents extension of reads beyond certain threshold and at the same time have an unlimited potential number of antibodies in the sample. (B) Schematics of experiment demonstrating the ‘‘activator’’ method and its robustness. Each antigen of a set of 22 surface markers is redundantly detected by three CODEX tag conjugates of the same antibody. The first conjugate is detected during panel 1 rendering, second – during panel 2 etc.. Thus the signal for same st th th antigen is detected at different cycles (e.g., 1 ,13 , and 24 ). (C) Montage of a fragment of imaging field of the 36 cycles of CODEX used to render a mixture of 18 barcoded spleens (similar to design in Figure 2A). Cycles N,N+12 and N+24 all three of which render same pair of antigens are shown per tile for all 11 pairs of antigens (see annotation in the black rectangle of each tile). Figure S3. Types of Samples in MRL/lpr Dataset, Related to Figure 5 MRL/lpr dataset has 9 samples: 3 control wild-type BALBc spleens (BALBc 1,-2,-3 and 6 MRL spleens MRL 4,-5,-6,-7,-8,-9). Based on disintegration of marginal zone as measured by frequency of marginal zone macrophages (MZM’s, – see black asterisk on Figure 5B and yellow arrow in this figure pointing to the area where CD169 positive (red) rim of MZMs is expected to be observed) and accumulation of double negative T cells expressing B220 B cell marker (B220 DN T cells – see red asterisk on Figure 5B) MRL spleens were grouped into early (MRL 4,-5,-6), intermediate (MRL 7,-8), and late (MRL 9) types. Early stage was represented by 3 MZM positive DN T cell-low spleens. Two spleens represented the intermediate stage: MZM low DN T cell-low spleen (Int1) and MZM positive DN T cell-positive spleen (Int2). Late stage was represented by single MZM positive DN T cell-positive spleen. A single representative spleen is shown for each stage together with interaction matrix. Color represents odd ratios (observed frequency of interaction/ expected frequency of interaction). Figure S4. CODEX Pinpoints Splenic Location of Unique Cell Types, Related to Figures 3–5 (A) Distribution of CD4(+)MHCII(+) cells (marked with white circles) in BALBc #2 spleen stained with IgD (green) and CD90 (red) to indicate positions of B and T cells accordingly. (B) CD4 and MHCII expression in isolated mouse splenocytes gated negative for all CODEX panel markers and in addition 120 g8 (lineage depletion with BD 558451 and dump channel for FITC conjugated or biotinilated antibodies corresponding to the antigens stained with CODEX panel were used for negative gating) except CD4, MHCII, CD45 and CD44. (C) CD4(+)MHCII(+) cells within the gate shown in (B) were sorted out and subjected to microarray analysis. CD4 T cells, CD8 T cells, bulk B cells and Conventional CD11c positive dendritic cells were co–sorted as a control. Expression of Lti signature genes (two individual signature sets as inferred in (Robinette et al., 2015)) in sorted cells. (D and E) CD11c+ B cells (age associated B cells (ABCs) in normal nd M/lpr spleens. ABCs have been shown to be a key participant in the triggering of certain autoimmune responses (Rubtsova et al., 2017, Rubtsov et al., 2011)) their splenic location has not been previously described in the literature. We observed ABCs to tightly associate with conventional dendritic cells (cDC) and occupy a distinct peri-follicular space in the boundary between PALS and B-zone. Interestingly, these cells diminished in numbers and redistributed toward intra-follicular space in the MRL/lpr spleens. (legend continued on next page) (F and G) Co-distribution of B220 and TCRb in isolated splenocytes of normal (BALBc) and autoimmune (MRL/lpr) mice. Gate in (G) points to significant (13%) presence of B220+ DN T cells in MRL spleen. (H–J) Thread like arrangement of CD8 T cells (purple, annotated with V-letter) has been noticed in PALS of splenic samples across dataset. To examine potential mechanisms driving these structures CD8 Tells and B220 positive B cells were sorted individually from BALBc spleen (I) and later combined in flat bottom microwell plates and mixed at 37C in culture medium. After mixing cells were stained for B220 (green) and CD8a (red) and imaged (J). Thread like structures similar to what was observed in spleen were detected. (K) Heatmap showing average frequencies of cell types (rows of heatmap) in the ring of index cell neighbors (see schematics on the right) for all niche clusters (0-99 in columns). (L) Heatmap shows how different cell types (in rows) are distributed between niches (in columns). Figure S5. ‘‘Cell Passports’’ of Selected Cell Types Identified in Normal and MRL Spleens, Related to Figures 3F and 5B (A) Diagram of per cycle markers for CODEX cycle montages in B,C and D. (B, E, and H) High resolution montage of CODEX cycles with cells of interest (CD11c(+) B cells) marked with yellow crosses is shown in (B). Low resolution montage of distribution of cells of interest (marked with white circles) in all imaged samples is shown in (E). Average expression profile of all markers in the cells of the selected cell type is shown in (H). (C, F, and I) Same for CD4(+)MHCII(+) cells. (D, G, and J) Same for CD106(+)CD16/32(-)Ly6C(+)CD31(+) cells. More examples of ‘‘cell passports’’ can be found in associated online repository (see STAR Methods). Figure S6. Cross-Tissue and Cross-Samples Distribution of Interacting Cell Pairs for Selected Types of Cell-to-Cell Interactions, Related to Figure 5C Interacting cell pairs are marked with white and cyan circles on the montage of IgD CD90 (B cell and T cell markers) staining of every sample of the dataset. Due to cell proximity in most cases cyan circles practically completely overlay white. (A and B) Interaction of CD4 and CD8 T cells with ERTR stroma (change in odds ratio score correlates with change in interaction count). (C–E) Interaction of granulocytes with CD4 T cells, dendritic cells and erythroblasts. (F and G) Interaction of erythroblasts with stromal and B220(+) DN T cells. Interactions in (C-G) scored as increased in early MRL/lpr (4,-5,-6) as compared to BALBc spleens (FDR of t test on normalized interaction counts between conditions < 0.05, difference in interaction counts > 0). (H and I) Interactions of B220(+) DN T cells with CD8 T cells and stromal cells. These interactions scored as increased in intermediate and late MRL/lpr (7,-8,-9) as compared to early MRL/lpr spleens (FDR of t test on normalized interaction counts between conditions < 0.05, difference in interaction counts > 0). More examples of cell type pairs with change in interactions across dataset can be found in associated online repository (see STAR Methods). Figure S7. The Fluidic Setup and Stage for Running CODEX Experiments, Related to STAR Methods (A) General diagram of robotic fluidic setup used in this study. CODEX experiments are done in an open flow cell, which can be imaged in any inverted mi- croscope. Six solutions have to be programmatically delivered and removed from the flow cell, which in the meantime sits in spatially defined position in the imaging system. A combination of 6-channel Tecan syringe pump equipped with 250ul syringes and USB-relay driven vacuum valve was used for iterative solution delivery and removal. Imaging was performed in Keyence BZ-X710 fluorescent microscope configured with 3 fluorescent channels (FITC, Cy3, Cy5) and equipped with Nikon PlanFluor 40x NA 1.3 oil immersion lens. Insets show photographs of actual microscope stage and fluidics robot. (B) Detailed 3D model of CODEX stage used in experiments. A metal insert was machined to be compatible with either ASI (Advanced Scientific Instrumentation) or Keyence 3d stages. Disposable (one per experiment) acrylic platform with a circular cutout in the middle was custom designed and lasercut such that it could be attached to the metal stage insert. Before multicycle run the coverslip with a sample was glued to the acrylic base which produced an open flow cell. As opposed to closed, open flow cell design ensures efficient (99.9%) and rapid solution exchange that is critical for CODEX protocol. (C) An exemplary photograph of full CODEX setup when attached to an inverted confocal microscope. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Cell Unpaywall

Deep Profiling of Mouse Splenic Architecture with CODEX Multiplexed Imaging

Loading next page...
 
/lp/unpaywall/deep-profiling-of-mouse-splenic-architecture-with-codex-multiplexed-j80eV0AV0q

References

References for this paper are not available at this time. We will be adding them shortly, thank you for your patience.

Publisher
Unpaywall
ISSN
0092-8674
DOI
10.1016/j.cell.2018.07.010
Publisher site
See Article on Publisher Site

Abstract

Resource Deep Profiling of Mouse Splenic Architecture with CODEX Multiplexed Imaging Graphical Abstract Authors Yury Goltsev, Nikolay Samusik, Julia Kennedy-Darling, ..., Gustavo Vazquez, Sarah Black, Garry P. Nolan Correspondence gnolan@stanford.edu In Brief A DNA barcoding-based imaging technique uses multiplexed tissue antigen staining to enable the characterization of cell types and dynamics in a model of autoimmune disease. Highlights d Autoimmunity analyzed by multiplexed DNA-tagged antibody staining (CODEX) d CODEX data reveal pairwise interactions and niches changing with disease d First tier of neighbors significantly impacts marker expression in the index cells d Changes in splenic morphology correlate with shifts in cell frequencies Goltsev et al., 2018, Cell 174, 968–981 August 9, 2018 ª 2018 The Author(s). Published by Elsevier Inc. https://doi.org/10.1016/j.cell.2018.07.010 Resource Deep Profiling of Mouse Splenic Architecture with CODEX Multiplexed Imaging 1,2 1,2 1 1 1,3 1 Yury Goltsev, Nikolay Samusik, Julia Kennedy-Darling, Salil Bhate, Matthew Hale, Gustavo Vazquez, 1 1,4, Sarah Black, and Garry P. Nolan * Department of Microbiology and Immunology, Baxter Laboratory in Stem Cell Biology, Stanford University School of Medicine, Stanford, CA, USA These authors contributed equally Present address: Smart Tube, Inc., 922 Industrial Avenue, Palo Alto, CA 94303, USA Lead Contact *Correspondence: gnolan@stanford.edu https://doi.org/10.1016/j.cell.2018.07.010 SUMMARY tation and emission spectra makes it hard to image more than 4–5 fluorophores with conventional fluorescent microscopy. A highly multiplexed cytometric imaging approach, Yet considerably more surface markers are needed for precise identification of cellular subsets and their activation state (Chat- termed co-detection by indexing (CODEX), is used topadhyay and Roederer, 2012). Approaches have been devel- here to create multiplexed datasets of normal and oped to overcome such limitations (Schubert et al., 2006; Gerdes lupus (MRL/lpr) murine spleens. CODEX iteratively et al., 2013), but these protocols have required multiple stain/ visualizes antibody binding events using DNA barc- strip/wash cycles of the antibodies that can be time consuming odes, fluorescent dNTP analogs, and an in situ or lead to sample degradation over the iterations. polymerization-based indexing procedure. An algo- The technique described here (CODEX, for CO-Detection by rithmic pipeline for single-cell antigen quantification indEXing) extends deep phenotyping capabilities of flow and in tightly packed tissues was developed and used mass cytometry (Spitzer et al., 2015; Bendall et al., 2011)to to overlay well-known morphological features with most standard three-color fluorescence microscope platforms de novo characterization of lymphoid tissue architec- for imaging of solid tissues. Accurate highly multiplexed single- ture at a single-cell and cellular neighborhood levels. cell quantification of membrane protein expression in densely packed lymphoid tissue images (which was once deemed We observed an unexpected, profound impact of the impossible [Gerner et al., 2012]) was achieved using polymer- cellular neighborhood on the expression of protein ase-driven incorporation of dye-labeled nucleotides into the receptors on immune cells. By comparing normal DNA tag of oligonucleotide-conjugated antibodies, combined murine spleen to spleens from animals with sys- with an image-based expression estimation algorithm. Auto- temic autoimmune disease (MRL/lpr), extensive matic delineation of cell types from multidimensional marker and previously uncharacterized splenic cell-interac- expression and positional data generated by CODEX enabled tion dynamics in the healthy versus diseased state deep characterization of cellular niches and their dynamics dur- was observed. The fidelity of multiplexed spatial cy- ing autoimmune disease both for major and rare cell types popu- tometry demonstrated here allows for quantitative lating mouse spleen. A rich source of multivariate data has been systemic characterization of tissue architecture in generated and provided for the community to further efforts in normal and clinically aberrant samples. developing approaches for image analysis, tissue architecture mapping, and rare cell-type detection. INTRODUCTION RESULTS Dramatic immune tissue re-organization has been seen in lupus erythematosus, where a variety of organs (from skin, to kidney, Single Base Primer Extension Enables Multiplexed and other body organs) can be targeted in relapsing-remitting Antigen Staining flares. One example of such reorganization is pronounced DNA provides an ideal substrate for designing molecular tags lymphadenopathy and splenomegaly observed in lupus models due to its combinatorial polymer nature. An indexable tagging (Lieberum and Hartmann, 1988; Jacobson et al., 1995). Using system whereby tags are iteratively revealed in situ by a stepwise mice with MRL/lpr genotype (Kanauchi et al., 1991), we sought enumeration procedure was designed. Antibodies (or other to systematically characterize microenvironment and cell inter- affinity-based probes) are labeled with uniquely designed oligo- actions associated with changes in immune organ architecture nucleotide duplexes with 5 overhangs that enable iterative step- and the progression of autoimmune disease. To this end, we wise visualization (Figure 1A; Video S1, part 1). Cells are stained devised a multiplexed microscopy technique that allows a pre- with a mixture of all tagged antibodies at once. At each rendering cise mapping of cell types in tissues. Significant overlap in exci- cycle, the cells are exposed to a nucleotide mix that contains one 968 Cell 174, 968–981, August 9, 2018 ª 2018 The Author(s). Published by Elsevier Inc. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Figure 1. Sequential Primer Extension on Samples Stained with DNA-Barcoded Antibodies Enables Unlimited Level of Multiplexing (A) CODEX schematic diagram. (B and C) Mouse spleen cells were fixed and co-stained with conventional TCR-b Ax488 antibody and CD4 antibody conjugated to CODEX oligonucleotide duplex as in first round of (A). After staining, cells were incubated in extension buffer with dG and dUTP-Cy5 either without (B) or with (C) Klenow exo- polymerase. Note that TCR-b-positive T cells in (B) and (C) are indicated by Ax-488 staining. Dependent upon the addition of Klenow, TCR-b-positive CD4 positive T cells are seen as a Cy5 positive subset of TCR-b-positive T cells in (C). (D) Spleen cryosection stained with B cell-specific B220-APC (red) and T cell-specific TCR-fluorescein isothiocyanate (FITC) (green) show mutually exclusive staining pattern in the marginal area between B cell follicle and the white pulp. (E) Spleen cryosection stained with CODEX DNA-tagged B220 (red) and CODEX DNA-tagged TCR-b (green) shows staining similar to the one observed with regular antibodies in (D). (F) Spleen sections were co-stained with regular B220-FITC and two antibodies (ERTR7 and CD169) tagged with cycle 1 CODEX DNA duplexes. Localization of marginal zone CD169 positive macrophages in the area between the ERTR7 positive splenic conduit of the white pulp and the B220 positive follicular B cells (D) as reported previously has been observed. See also Figure S1 and Video S1, part 1. of two non-fluorescent ‘‘index’’ nucleotides and two fluorescent is paused at the indexing position by omitting one of the indexing labeling nucleotides. The index nucleotides fills in the first index (walking) bases from the labeling mix (as done in this study) or position across all antibodies bound to the cells. However, the potentially by use of reversible terminators (Video S1, part 2). DNA tags are designed such that only the first two antibodies Importantly, the system enables multiplexed tissue imaging are capable of being labeled with one of the two fluorescent analysis by means of a standard fluorescence microscope. dNTPs—and only if the index nucleotide was previously incorpo- To test the premise of the system, isolated mouse spleen cells rated. Those two antibodies are then imaged by standard fluo- were incubated with a CD4 antibody conjugated to an indexing rescence microscopy. Then the fluorophores are cleaved and oligonucleotide duplex (as represented by Ab1 in Figure 1A). In washed away, and the sample is ready for the next cycle where this trial experiment, T cell receptor b (TCR-b)-Alexa 488 was a different indexing nucleotide is used. At the end of the multi- used as a counterstain. A single round of primer extension was cycle rendering protocol each pair of antibodies is visualized at done with a mix of unlabeled dGTP and dUTP-ss-Cy5. A cell a known, pre-defined cycle of the indexing protocol, and the population positive for both CD4 and TCR-b was observed by multiparameter image can be reconstructed. The polymerase flow cytometry. Observation of this population was dependent Cell 174, 968–981, August 9, 2018 969 on the addition of Klenow DNA polymerase to the reaction (Figure 2). Indeed, use of this compensation method resulted mixture (Figures 1B and 1C) indicating specific antibody in a considerable (approximately 2-fold) reduction of spillover rendering by primer extension. Similarly, in tissue sections, signal (especially pronounced for CD4/CD8a co-distribution— CODEX tag-conjugated antibodies produced lineage-specific Figures 1L and 2E). staining comparable to regular fluorescent antibodies (compare A 30-antibody panel was therefore designed to identify staining patterns of B220-CODEX and B220-APC in mouse splenic-resident cell types (lymphocytes, macrophages, micro- spleen, Figures 1D–1F). vessels, conduit system, and splenic stroma; Figure 3A; Table A simulated multicellular mix was produced by combining 30 S1) and applied to the cryosections of spleens from wild-type aliquots of mouse splenocytes barcoded with pan-leukocytic (3 spleens) and MRL/lpr mice (6 spleens) (Figure S3). Four major CD45 antibody labeled with one of 30 distinct CODEX tags (Fig- classic splenic compartments: red pulp, B cell follicle, PALS ure S1A). The visualization of the CODEX 15-cycle staining (periarteriolar lymphoid sheath), and marginal zone (MZ) (Fig- pattern showed even cycle-specific signals, with low back- ure 3B) could be easily discerned in CODEX imaging data (Fig- ground (average signal to noise 85:1), efficient (98%) release ure 3A). A total of 734,101 segmented cells were identified, of fluorophore by inter-cycle TCEP (Tris(2-carboxyethyl)phos- and by means of X-shift clustering (see STAR Methods) their phine hydrochloride) cleavage and no signal carryover between expression profiles were grouped into 27 broadly defined pheno- cycles (Figures S1B–S1D and S1F). Linear regression analysis typic groups (Figures 3C–3F; Table S2) most of them matching revealed low signal deterioration (at 0.79% per cycle) and to known cell types. Compared to CyTOF data on splenocytes acceptable background (starting at 1.1% and increasing at isolated from non-enzymatically homogenized spleen, CODEX 0.06% per cycle, Figures S1E and S1F). in situ analysis produced a similar distribution of cell counts for To further reduce the signal loss associated with accumulation major cell type. Yet being a non-disruptive technique CODEX of polymerization errors and to allow larger panels without identified larger numbers of resident and stromal cell types increasing the length of tagging oligonucleotides, an approach such as erythroblasts and F4/80 macrophages than CyTOF did based on primer-dependent subpanels was devised (Fig- (Figure 3E). Notably even rare computationally derived cellular hi – hi + ure S2A). The feasibility of this design expansion was tested by phenotypes (e.g., CD4 /CD3 /MHCII cells and CD11c B cells) staining mouse splenocytes with a 22-plex set of antibodies. closely matched the cell types previously observed in murine Each of the antibodies was conjugated to the three versions of spleen (LTi cells [Robinette et al., 2015] Figures S4A–S4C, the CODEX duplex tag—with same terminated top oligonucleo- S5C, S5F, and S5I and age-associated B cells [ABCs], Figures tide and three kinds of the tagging oligo (Figure S2B). Thus, every S4D, S4E, S5B, S5E, and S5H). antigen was detected thrice (bringing the overall number of detections to 66) and only after annealing of a panel-specific Pairwise and Combinatorial (i-niches) Statistics of Cell- activator oligonucleotide. We found that the signal for same to-Cell Contacts in Mouse Spleen antibody was consistent across the three primer-dependent To provide a high-level view of the cell-type interaction land- batches (Figure S2C). Thus, panel-activator design extends scape, the total counts of contacts between every pair of cell CODEX to a theoretically unlimited multiplexing capacity, types in the Delaunay neighborhood graph (Gabriel and Sokal, bounded only by the speed and resolution of the imaging pro- 1969)(Figure 4A and associated Mendeley dataset) for each cess itself and the time required for each imaging cycle. condition was determined. The specificity of cell-to-cell interac- tion was estimated from the ‘‘log odds ratio’’ metric (log-ratio of Benchmarking CODEX observed and expected probabilities of contacts between 2 cell To validate the quantitative performance of CODEX, cells freshly types) (Table S3). When visualized as heatmaps, this metric re- isolated from mouse spleens were co-analyzed by mass-cytom- vealed a significant non-random distribution of cells in the etry (CyTOF) and CODEX using identical 24-antibody panels spleen. In the majority of cases, cell types were either selectively (Table S1). Use of the same antibody clones and the same sple- associating or avoiding each other (red or blue on the heatmap) nocyte preparation ensured the validity of comparisons. Antigen pointing to prevalence of specific cell-to-cell interactions in co-expression signals from CODEX, obtained from image seg- shaping the spleen architecture. The major splenic anatomic mentation (see STAR Methods), were consistently similar to compartments were reflected in two large mutually exclusive CyTOF (Figure 2B, for direct comparability, both CODEX and clusters of positive associations, which appeared to correspond CyTOF data are plotted on a linear scale). to red pulp and the white pulp, respectively (indicated with black Consecutively, CODEX was applied to tissue sections. In rectangular outlines on Figure 3G). For example, a significant contrast to dissociated cells spreads (Figure 2A), cells in tissue positive association was observed between F4/80 macro- sections are adjacent to each other—with large fractions of phages and erythroid cells, as these cell types are both found membranes in direct contact (Figure 2C). Therefore, neighboring in the red pulp and are closely associated in so-called erythro- cells can contaminate each other’s signals during the quantifica- blast islands (Socolovsky et al., 2007). An avoidance of interac- tion phase (Figure 2D). To address this latter challenge, a novel tion was observed between T and B cells, reflecting concentra- linear algorithm for 3D positional spillover compensation was tion of these cell types in B cell follicles and PALS, respectively created. This algorithm is based on the same principles used (Figure 3G). in fluorescent spillover compensation in traditional flow cytome- Unexpectedly, a consistently high degree of association was try, except that our algorithm performs compensation between observed between the cells of the same phenotypic class physically adjacent cell based on their surface contact ratios (Figure 3G, red diagonal), suggesting that homotypic adhesion 970 Cell 174, 968–981, August 9, 2018 Figure 2. Accuracy of Surface-Marker Quantitation by CODEX (A) Microscopic image of mouse splenocytes stained with a 24-color antibody panel, showing one cycle of CODEX antibody rendering. Cell contours show the outlines produced by the cell segmentation algorithm (B) Comparison of single-cell expression data derived from dissociated mouse splenocytes on an identical 24-color panel using CODEX and CyTOF. (C) Example segmentation in a mouse spleen section based on combining nuclear and membrane (CD45) channel. (D) Graphical explanation of the algorithm for compensating the spillover between neighboring cells using a cell-by-cell compensation matrix. (E) Biaxial plots of segmented CODEX data acquired in mouse (BALBc) spleen sections. The presence of double-positive cells in the upper quadrant is used as an estimate of lateral signal bleeding explained schematically in (D). Three combinations of mutually exclusive lineage markers are shown to demonstrate the range of effect of the compensation algorithm on reduction of lateral signal bleed. See also Figures S4 and S5. constitutes a major force driving the architecture of immune tis- previously. We define here an indexed ‘‘niche’’ (i-niche) as a ring sue. This observation held true both for the major constituents of of cells (excluding the central, or here defined as ‘‘index’’ cell) in white pulp, T and B cells, as well as for rare cell types such as no specific circumferential order that are Delaunay neighbors of natural killer (NK) cells. Interestingly, even though CD8 and the index cell (Figure 4A and STAR Methods). We identified 100 CD4 T cells tended to mix in the PALS, their mutual distribution of the major i-niches (by K-means clustering) according to the was nonrandom and consisted of intertwined threads of homo- relative frequency of the identified cell types present in the ring typic cells (Figure S4H). Interestingly, as an aside, similar struc- of cells surrounding the index cell (Figures 4A, 4B, and 4F). tures could be reproduced in vitro by incubating heterotypic Most i-niches could be readily mapped into one of major mixtures of sorted splenic cell populations (Figures S4I and anatomical compartments of the spleen (B cell follicle, PALS, S4J). These data suggest that homotypic cell association marginal zone, or red pulp, per Figure 4G). In most cases, any might be an important driver of the white pulp substructure given i-niche resided within a single anatomical compartment and is worth further investigation. (although several i-niches were observed in more than one The precision in situ cytometry analysis of CODEX data al- compartment), and every splenic compartment was populated lowed enumeration of cellular contexts in a manner not possible by many i-niches (Figure 4G). Cell 174, 968–981, August 9, 2018 971 Figure 3. CODEX Analysis of Mouse Spleen Cryosections Co-stained for 28 Antigens (A) Three collated images on the left correspond to the legend of antibody renderings per cycle, gross morphology photograph of MRL/lpr (left), and normal (right) spleen embedded in optimal cutting temperature (OCT) block prior to sectioning. Green corresponds to antibodies rendered by extension with dUTP-Cy5; red, dCTP-Cy3 On the right, collage of the CODEX multicycle data for normal spleen (BALBc-2) and early MRL/lpr spleen (MRL/lpr 4). All images are derived from a single scan with a 403 oil objective of an area covered by 63 tiled fields. (B) Schematic diagram of major known splenic anatomical subdivisions drawn based on cell distribution in BALBc-1 replicate. (C) An exemplary profile of Vortex cluster (B cells) used for manual matching of clusters to known cell types. (D) Minimal spanning tree (MSP) built for all clusters identified by Vortex analysis. On the left middle and right panels, the MSPs are colored by expression levels of B220, TCR, and CD71 accordingly to indicate location of B cells T cells and erythroblasts on the tree. (E) Circle chart showing for several major cell types their fraction of total cells as identified by CODEX analysis of splenic tissue and CYTOF analysis of isolated BALBc splenocytes (F) Post-segmentation-derived diagram of identified objects (cells) colored according to cell types in BALBc-1 replicate. Full-size diagrams for every tissue analyzed in this study are available online (see STAR Methods) (G) Average cell-type to cell-type interaction strength heatmap for BALBc samples. Color from blue (<0) to white (around 0) to red (>0) indicates log of odds ratio of interaction (ratio of observed frequency versus expected frequency of interaction). The rows and columns are in the same order (annotation on the right). Black outlines indicate two largely exclusive mega-clusters of cross-interacting cell types loosely matching the cell types populating the red and the white pulp. See also Figure S5 and Video S2. 972 Cell 174, 968–981, August 9, 2018 Figure 4. Unbiased Identification of i-niches in Multidimensional CODEX Data (A) On the left diagram explaining the terminology used for defining i-niche (a ring of first tier neighbors for central cell). On the right Delaunay triangulation graph used for identification of first tier of neighbors for every cell. (B) Heatmap depicting frequency of cell types in 100 types of i-niches identified by K-means (K = 100) clustering of all index cells in the dataset (each cell is an index cell for its i-niche) based on frequency of different cell types in the first tier of neighbors. The color indicates the average fraction of corresponding cell type in the the i-niche. (C) An example of marginal zone and follicular (B-zone) B cells defined by residence in distinct i-niches (e.g., marginal zone i-niche includes a marginal zone macrophage marked by letter H and green color). Positions of B cells in each i-niche is marked with red circles over the schematic of BALBc spleen. (D and E) Two heatmaps from top to bottom show average expression of selected surface markers measured in a central cell across 100 i-niches (same left to right order as in B) when central cell is B cells (D) or CD4 T cell (E) accordingly. The color indicates the relative level of surface-marker expression as measured across dataset. Gray columns indicate absence of cells in corresponding niches. Two orange rectangles over top heatmap indicates position of i-niches with high CD35 (containing FDCs and marginal zone macrophages). Cyan rectangle shows location of family of i-niches with high content of F4/80 macrophages and low B220 and CD19 in central B cells. Purple rectangle indicates family of i-niches enriched with ERTR-7 positive stroma. Below top heatmap, location of selected i-niches shown in (E) are indicated. Over bottom heatmap, yellow rectangle indicates the family of i-niches with dominating presence of B cells. Two green rectangles indicate family of niches with high levels of CD90 and CD27 in the index CD4 T cells. (F and G) Abundance of 100 i-niches in normal spleen (top bar graph) (F) and relative distribution of i-niches (G) between splenic histological subdivisions (PALS, red pulp, marginal zone, and B-zone) shown as a heatmap. To illustrate a variety of tissue distribution pattern by i-niches an overlay of selected i-niches over a schematic of normal spleen (BALBc-1) is shown. Heatmap color indicates fraction of corresponding i-niche per splenic anatomic subdivision. (H) Top right shows a biaxial plot of flow data for CD79b and B220 measured in isolated splenocytes. Top left shows levels of CD79b and B220 in central B cells as measured across all 100 i-niches. To illustrate i-niche-dependent variability of surface-marker expression, images of central cells (marked with red cross) with levels of surface marker indicated in pseudocolor palette are shown for selected exemplary i-niches in the bottom panels. See also Figures S4K and S4L. Cell 174, 968–981, August 9, 2018 973 In our definition, the index cell in the center can be any cell. levels on splenic B cells could be, to a large degree, accounted Thus, i-niches enabled subsetting the common cell types based for by the niche composition around those B cells—and that on cellular context (microenvironment). For example, B cells sur- the expression levels on these cells might be influenced by (or rounded by only B cells (red arrow Figure 4B, i-niche #96) can be influences) the cells in their immediate surrounding. seen primarily in the follicular zone B cell region (Figure 4C, left The overall utility of the i-niche in determining any given sur- panel), while presence of CD169 positive marginal zone macro- face-marker expression value for an index cell was evaluated phages mapped the B cells in such i-niches to the marginal zone by constructing a linear regression model of marker expression (red arrow Figure 4B, i-niche #33, Figure 4C, right panel). In the using both the cell-type identity and the i-niche constituency in case of T cells, CODEX data enabled precise selection of an a two-featured variable model (the other variable being the cell- important migratory subset of T cells known to be residing in type identity). Notably, adding the i-niche information as a depen- ERTR7 enriched niches (Burrell et al., 2015) (in Figures 4B and dent variable significantly improved the fitness of the model for all 4E see the purple rectangle indicating a family of niches where markers (Table S4) with highest improvement F-values for CD90, index T cells contact ERTR7 stroma; as well as Figures S6A B220, CD21/35, and ERTR7 and the lowest prediction rates for and S6B). Taken together, we see that surface-marker expres- Ly6G, CD5, CD11b, CD5, and TCR-b. Thus, the high variability sion alone is insufficient to associate many cell subsets with a in B220 expression levels are highly related to the i-niche in which given tissue subcompartment (e.g., CD4 T cells can be found the B cell resides. In other words, B220 expression levels can be both in the PALS and in the red pulp). However, i-niche designa- location specific and are dependent on the i-niche partners. As a tion does provide such mapping data (most of i-niches were en- counter point, the data also show that i-niche does not reliably riched within a specific splenic subdivision Figure 4G) and as predict expression of other proteins, such as CD5 or TCR-b, such T cells associated with ERTR7 positive stroma in fact the expression levels of these receptors is relatively constant localize primarily to the red pulp (Figures S6A and S6B). This rai- across the i-niches (Figure 4E). This result quantitatively demon- ses interesting questions—can new cell types, or functional sub- strates that the i-niche (neighbors) determines a significant pro- sets, be discerned by this approach? What is the frequency of a portion of variance in the expression of certain markers. Overall, repeated i-niche structure that must be observed to suggest a we observed that that many splenic cell types populate a wide va- function? And what would constitute a proof that a given i-niche riety of i-niches (Figures S4K and S4L), suggestive of a multiplicity corresponds to a new cell type or functionality? of functional state for any given immune cell type. Further, tissue locale (i-niches) is a powerful indicator of potential differential Cell-Surface-Marker Expression Depends on Local function (to the extent tissue locale drives function), and these Neighborhood deterministic changes in surface marker protein expression are One approach to address these latter questions is to consider surrogate indicators of this locale or function. the phenotypes of the index cell in various i-niches. We observed that for several index cell types (e.g., for B and T cells) there was Changes in Splenic Composition Associated with significant biasing of the surface-marker expression depending Disease Progression on the i-niche in which the index cell resides (e.g., see selected A comparable region of spleen was visualized by CODEX for 3 cases marked with colored rectangles above the heatmaps in normal BALBc spleens, and 6 spleens from MRL/lpr mice. Image Figures 4D and 4E). To assure that it was not a quantitation arti- segmentation revealed strong variation in cell counts between fact, we mapped back to the image the index B cells where the the norm and the disease (Figure 5B) for most (19 out of 27) of levels of CD79b (a co-activator chain of the B cell receptor com- the cell types identified by X-shift clustering. Examples include plex) and B220 (a splice isoform of CD45 membrane phospha- a dramatic increase in CD71 erythroblasts (green cells on Fig- tase) would be niche dependent. We found that the index B cells ure 5A maps), a reduction in numbers of B cells and FDCs (follic- int lo + that were B220 , CD79b (i-niche ‘‘59’’) resided on the ular dendritic cells), and increases in so-called B220 DN T cells boundary between the PALS and the follicles (Figure 4H, image (CD4/CD8 double-negative B220 T cells), which have been pre- lo montage on the bottom). Index B cells that were B220 , viously characterized as a hallmark of the MRL/lpr progression int CD79b (i-niche ‘‘91’’) were mostly found in the red pulp. And, (Koh et al., 1995) which could also be identified by fluores- int/hi hi B cells that were B220 , CD79b (i-niche ‘‘76’’) were yet cence-activated cell sorting (FACS) (Figures S1F and S1G), different again and were found at the boundary of the red pulp thus ruling out the possibility that this unusual cell type being a and the follicles. When measured in cell suspensions, CD79b result of image segmentation errors. These and other changes is co-expressed with B220 at various ratios (see CyTOF plot sin- were used to broadly classify the MRL/lpr spleens into early, in- gle-cell splenocytes Figure 4H, top-right panel). Such a distribu- termediate, and late disease stages (Figure S3). tion of expression is sometimes attributed to staining variability, Despite consistent presence of ‘‘homotypic interactions’’ di- measurement noise, or a simple lack of understanding of the agonal and larger cell-adjacency clusters corresponding to red underlying biology. However, as seen on Figure 4H, upper-left and white pulp in odds ratio heatmaps across disease (Fig- panel, there is a non-random pattern of CD79b and B220 expres- ure S3), a deeper statistical analysis revealed many disease- sion across the central cell of the corresponding i-niches, and, associated changes in frequency of contacts between cell types depending on the B220/CD79b levels, the i-niches (the central (see Table S3). Among the changes we observed an increase in – + cells) map to specific regions in the splenic architecture (Fig- interaction between B cells and CD4 /CD8 conventional den- ure 4H, lower panel). These observations suggest that the dritic cells (cDCs) in the early MRL/lpr spleen compared to spread of the CD79b-B220 levels as well as of other marker normal (Figure 5C, left panel), suggesting an increase in B cell 974 Cell 174, 968–981, August 9, 2018 Figure 5. Autoimmune Disease Drives Changes in Splenic Composition and Cell-to-Cell Interactions (A) Post-segmentation diagrams of all objects (cells) colored according to cell types (see color map in Figure 3F) for all normal and MRL/lpr tissue sections imaged in the study. Full-size diagrams are available for every tissue analyzed in this study are available online (see STAR Methods). (B) Stacked bar graphs show dynamics of cell counts across dataset for manually annotated Vortex clusters (cell types on the left) across progression from normal to afflicted spleen. Colored bar sections indicate fraction of the total cells as detected at a particular stage/samples (1–9 annotation on the top). Cell types were split into four types according to the dynamics of counts across dataset as represented by average relative (normalized to 1) count; see line graphs on the right; x axis corresponds to stage/sample id. (C) Two examples of change in cell-to-cell interaction frequency during disease progression between the B cells and dendritic cells in normal and early MRL/lpr spleen and between B220 DN T cells and CD4 T cells during progression from early MRL/lpr to intermediate. (D) Co-distribution of odds ratio log fold [log(odds ratio in early MRL/lpr) – log(odds ratio in BALBc)] on x axis and change in counts of interactions for early MRL/lpr versus control (BALBc) comparisons (on y axis). (E) Co-distribution of cumulative cell-frequency change [celltype1 freq. change + celltype2 freq. change] on x axis and change in counts of interactions for early MRL/lpr versus control (BALBc) comparisons (on y axis). (F) Bar graph showing chi-square values across conditions computed for odds ratio and direct interaction counts. See also Video S2 and Figure S6. activation. We also observed a higher interaction frequency of S6G). In the intermediate- and late-stage MRL/lpr spleens, there granulocytes with T cells (Figure S6C), dendritic cells (Fig- was a significant increase in interaction of B220 DN T cells with + + ure S6D), and erythroblasts (Figure S6E), and a higher number CD4 T cells (Figure 5C, right panel), CD8 T cells, erythroblasts, of contacts between erythroblasts and various kinds of stromal and a variety of other cell types compared with numbers of these cells, as well as B220 DN T cells (Table S3; Figures S6F and interactions in the early MRL/lpr stage (Table S3 and Figures Cell 174, 968–981, August 9, 2018 975 S6G–S6I, see online resource, STAR Methods, for more exam- For the rest 24 of the 26 changing interactions mentioned ples). So, while there was no obvious gross rearrangement of above at least one of the cells of the pair was scored as signifi- the tissues, many homotypic and heterotypic cell-cell associa- cantly (FDR <0.05) changing the frequency across scored condi- tions were altered, prompting a key question: what are the tions (Table S3, last column of the ‘‘EarlyMRL vs BALBc control’’ main factors driving this disruption? spreadsheet). We therefore conclude that—at least in the diseased state of early-stage MRL/lpr—most of the change in Disease-Driven Change in Cell Counts Determines the counts of cell-cell interactions are driven simply by increases Frequency of Specific Cell-to-Cell Contacts or decreases in cell-type frequencies. What could be the drivers of changes in frequency of pairwise As an additional evidence, c statistics were used to compare cell-cell contacts? If the kinetics of pairwise cell contacts follows the total magnitude of changes in pairwise cell-type interaction a rate law, one possibility would be that modulation of specific matrices (total interaction count) versus changes in log-odds cell-to-cell interaction potential—or ‘‘attraction’’ (for which the ratio matrices (propensity for non-random interaction). The c odds ratio score was used as an estimate across this study)— deviation (sum of squares of Z-score-normalized values) was is the main driver. In other words, it would be expected that, computed for each disease matrix compared to the control. when the affinity of such an interaction goes up, the fraction In every case, the c values of cell-interaction matrices were of interacting cells of a given cell pair would increase. At larger than of the respective log odds ratio matrices of the the same time, even in the absence of change in cell-to-cell same biological sample (Figure 5F). This suggests that as affinity, the absolute number of the cell-cell pairs (defined here the cell-type frequencies change due to disease progression, the absolute numbers of interactions change dramatically as cell pair aggregates, or CPAs) and the number of interacting cell pairs should correlate with the frequencies of interacting whereas the frequency-normalized likelihoods of cell interac- cell types (analogous to concentrations in the rate law equation). tions change to a much smaller extent indicating a great degree Importantly, the latter scenario could be as biologically signifi- of robustness of the ‘‘design principles’’ of the splenic tissue and cant as the former. Finally, some of the cell-type contacts may that many of the more dramatic disease-associated variations be observed due to low cellular motility of randomly meeting occur primarily through the shift in cell numbers. cells. Such interactions would not produce spatially defined This analysis implies that the degeneration of the tissue integ- sub-splenic CPAs and would have and odds ratio close to 1. rity in the MRL disease largely follows dramatic changes in cell- The perturbation introduced to normal splenic composition type frequencies. At the same time, there were notable excep- with MRL/lpr genotype allowed us to examine the mechanisms tions from this trend, where the changes in observed cell-type implicated in transition from normal to diseased spleen. In short, pairing frequencies could be largely explained by shifts in the we found that, for most cell-cell pairs observed, the mutual cell-type interaction likelihoods (log-odds ratios). While further attraction (as quantified by the interaction log odds ratio) was work is required to determine which of these changes are instru- not the primary determinant driving the change in counts of inter- mental to the MRL disease state, our findings suggest that such acting cell pairs between the MRL/lpr and the norm. In Figure 5D, differential analysis can be applicable in other diseases and, we plot the change in counts of interactions of two cell types possibly, could be used to discover cell-type interactions that (e.g., A:B) between the MRL/lpr and the normal BALBc spleens. are targetable from a therapeutic standpoint. Each dot represents a pair of cell types. The value on the y axis is the difference in the total number of observed interactions be- Reorganization of Cells in Disease-Associated Tissue tween BALBc and MRL/lpr. The x axis shows the difference be- Substructures tween log odds ratios of interactions between the same condi- We cataloged the cell-cell interaction ‘‘connectivity’’ in a circular tions. There was no overall correlation observed (R = 0.058). correlation diagram. Rarely, if ever, there was any cell type found In contrast, we observed a correlation (R = 0.288) between adjacent to only one other type of cell. The highest degree of the cell count changes and the interaction changes (Figure 5E). connectivity was observed for the most abundant cell types In agreement with those observations, we saw that out of the such as B cells in normal spleen and erythroblasts (Figure 6A) 26 top-scoring (false discovery rate [FDR] <0.05 and change in early MRL/lpr. This high connectivity in turn led to large effect in absolute interaction counts >150) cell-type pairs of this on i-niches caused by changes in cell numbers associated with cross comparison, only 2 showed corresponding significant progression of disease from normal to autoimmunity. The most (FDR <0.05) change in odds ratio score. Curiously, these two in- dramatic changes in cell frequencies were the increase in eryth- teractions with a modest 1.5 times increase in interaction count roblasts in the early MRL/lpr and the emergence of B220 DN and, concomitantly, a 0.8 increase in log odds ratio score were T cells in late MRL/lpr—which were associated with the appear- the ones between the CD4 or CD8 T cells and ERTR7 stroma ance of novel i-niches relative to the normal spleen (spatial (see Table S3, rows 6 and 7, and Figures S6A and S6B). Visually localizations of B220 DN T cell-dominated i-niche 18, erythro- they appeared as persistent co-clustering of T cells with ERTR7 blast-driven i-niche 29, and B cell-rich i-niche 96 are shown on stroma despite the overall drop of T cell numbers in the ‘‘early’’ Figure 6B, and their cell-type composition is shown on heatmap MRL/lpr samples. Curiously, ERTR7 positive fibers of splenic on Figure 6C). A corollary to this is the question of whether the stroma as well as ERTR7 protein itself were recently shown to presence of these cells, and new i-niches dependent on these be critically involved in T cell trafficking (Burrell et al., 2015), sug- cells, somehow changed the observable biology of the cells gesting that this increase in the spatial association could be they contact? We found some examples supporting that, reflective of the T cell activation. whereby the proximity of CD4 T cells to B220 DN T leads to 976 Cell 174, 968–981, August 9, 2018 Figure 6. Differential Effect of Disease over i-niche Presence across Dataset (A) Cell-interaction networks built for BALBc early MRL/lpr and late MRL/lpr based on the number of contacts observed between two cell types (only connections with more then 150 interactions per sample are shown on the diagrams). Thickness of connection correlates with number of contacts size of the node indicates number of cells per condition. (B) Evolution of i-niche abundance across dataset. Selected three i-niches (marked above heatmap in C depicting i-niche composition) differentially represented across dataset (changing between norm and disease) are shown. Yellow circles overlaid over blank rectangles corresponding to imaged area indicate location of i-niche. (C) Top heatmap shows frequencies of B220 DN T cells, erythroblasts, and B cells in the i-niche rings. Line above top heatmap indicates the composition of i-niches 18, 29, and 96 described in (B). Color scheme is the same middle heatmap and indicates expression of selected markers when the i-niche central cell is an erythroblast—primarily to show that CD27 is not expressed on erythroblasts in the vicinity of B220 DN T cells. Bottom heatmap indicates expression of selected (legend continued on next page) Cell 174, 968–981, August 9, 2018 977 CD4 T cell activation in spleens of MRL/lpr mice: Figure 6C The neural network highlighted the regions in each multipa- shows increased levels of CD27 expression in CD4 T cells pre- rameter spleen image that corresponded to the disease state sent in i-niches dominated by B220 DN T cells (Figure 6C, red (Figure 7B), despite having seen no images from these spleens circle). during training. To investigate the specific features learned by Other cell types noticeably changed their characteristic distri- the neural network, the cell-type compositions of the regions bution and their propensity to engage, or evade, specific cell-to- identified as diseased versus those regions identified as normal cell contacts (as estimated by odds ratio score) during disease were compared. There was significant enrichment of several cell + – + progression. For example, cells of CD106 CD16/32 Ly6C types in these regions (Figure 7C). Although some cell types en- + + CD31 phenotype were randomly distributed in the red pulp riched in diseased regions, for example, B220 DN T cells, were of normal spleens but were found to aggregate in the areas present only in the diseased tissue, the most highly enriched cell + – proximal to the marginal zone of the MRL/lpr white pulp (Fig- type (CD4 /CD8 cDCs) was present in both the disease state ures S5D, S5G, and S5J). This re-distribution correlated with and the healthy state. erythroid proliferation and reduced odds ratio score for the inter- To assess the specific contextual changes recognized by the + – + + + – action of CD106 CD16/32 Ly6C CD31 and erythroblasts in neural network, the local neighborhoods of the CD4 /CD8 cDCs lupus spleens (Table S3). that the neural network found to be enriched in MRL/lpr regions were analyzed. In these neighborhoods, we observed a signifi- + – Automatic Definition of Disease-Associated Areas in cant enrichment of other CD4 /CD8 cDCs, as well as signif- + + – – Tissue Architecture icant depletion of CD106 /CD16/32 /Ly6C /CD31 stromal cells (FDR <10 ). This suggests that the neural network had identified As noted, the analysis reveals that the development of the auto- + – immune disease in mice (as exemplified by MRL/lpr lupus) is an altered context for CD4 /CD8 cDCs (distant from stromal associated with vast rearrangement of normal spleen architec- regions) as a key descriptor for the disease. Thus, the neural ture, which is likely to cause loss of cell-cell contexts normally network approach described here enabled both automatic clas- hosting the cells crucial for proper splenic function, as well as sification of samples according to disease state and an auto- the observed emergence of novel i-niches that are not found matic identification of high-dimensional regions of interest and in the normal BALBc spleen. Additionally, certain i-niches corresponding cellular niches. were sequestered to specific anatomic compartments of the spleen, which allowed us to use such i-niches as reference DISCUSSION points to quantitatively monitor high-order morphological changes. The i-niches that in normal spleen were localized to Here the feasibility of polymerase-driven highly multiplexed visu- one distinct compartment (more than 90% of central cells reside alization of antibody binding events to dissociated single cells as within a particular splenic compartment) were used to evaluate well as tissue sections (CODEX) was demonstrated and bench- the dynamics of splenic cells associated with progression of marked. Critically, CODEX enables co-staining of all antigens autoimmune disease (Figure 7A, middle heatmap). This analysis simultaneously with the staining iteratively revealed by primer confirmed the dissipation of the marginal zone starting from extension cycles wherein no diminution of epitope signal detec- early stages of MRL/lpr and revealed a progressive distortion tion was observed. A consistent performance of CODEX in co- of PALS. Curiously, depending on whether a i-niche was based detecting up to 66 antigens was demonstrated, and the ‘‘activa- on F4/80 macrophages or primarily contained erythroblasts, the tion primer’’-based extension of the system could enable a red pulp appeared to reorganize in the diseased tissue (Fig- potentially vast expansion of CODEX multiplexing capacity. For ure 7A, right heatmap), pointing to the fact that more than the current method, fresh-frozen tissue was used yet at a cost one compartment-specific niche is required to reliably trace of testing an extensively broader collection of clones we have the fate of specific anatomic compartments. In many cases, recently succeeded in adapting the procedure to formalin fixed the definition of subsets/morphological units constituting the paraffin embedded (FFPE) archival tissue (unpublished data). tissue is subjective, yet this study employed niches that were We believe this will open the large retrospective collection of algorithmically defined. Therefore, using niches as markers of FFPE samples from clinical cohorts to multidimensional cyto- morphology can quantitatively monitor the changes of high- metric analysis. order anatomic architecture. With the help of a simple automated fluidic setup (Figure S7), To automatically isolate the specific local combinations of the CODEX platform can be performed on most three-color fluo- expression patterns characteristic of the disease state, a fully rescence microscopes with a motorized stage and thereby convolutional neural network was trained to distinguish image enable conversion of regular fluorescence microscope into a patches from normal and MRL/lpr mice. The neural network tool for multidimensional tissue rendering and cell cytometry. operated by identifying, in each training image patch, the spe- CODEX completes a 30-antibody visualization in approximately cific areas that corresponded to the disease state. 3.5 hr. Modifications to the technology that increase the markers when the i-niche central cell is a CD4 T cell. The color schemes in these three heatmaps are the same as in heatmaps in Figures 4B, 4D, and 4E. Red oval outline pinpoints i-niches with elevated CD27. Note that these i-niches as indicated by top heatmap have B220 DN T cells as a prevailing component. Lower panels show examples of central cells in i-niches marked under the lower heatmap. i-niche 50 is an example of i-niche without B220 DN T cells. Central cell does not express high CD27. i-niches 42 and 44 have high frequency of B220 DN T cells and accordingly central cells express high CD27. See also Figures S4K and S4L. 978 Cell 174, 968–981, August 9, 2018 Figure 7. i-niches and Neural Nets Provide Unbiased Way for Disease Monitoring (A) Selected i-niches (green heatmap shows i-niche composition, color scheme same as in Figure 4B) were chosen based on high (>90%) presence per single histological subdivision (blue heatmap color scheme same as in Figure 4G). Abundance of these i-niches (brown heatmap, color indicates relative abundance of corresponding i-niche as measured across full dataset) was used to judge the preservation or decay of a histological splenic subdivision corresponding to selected i-niches. (B) Red color over blue rectangle indicates regions of interest (MRL/lpr-specific regions) predicted by neural network in entire spleen images. From top left, clockwise: BALBc #3, MRL/lpr #5, MRL/lpr #7, MRL/lpr #8. (C) Cell types enriched (FDR <0.1) in MRL/lpr-specific regions (in red in B) predicted by neural network. See also STAR Methods ‘‘Neural network training.’’ measurements per cycle, reduce the cycling time, faster imaging (MRL/lpr). Much like with conventional flow cytometry, CODEX methods such as light sheet microscopy, or an increased size of discerned all major cell types commonly observed in mouse the imaging the field of view offer potential opportunities for spleen. Moreover, application of X-shift phenotype-mapping al- increasing the depth and speed of the visualization process. gorithm (Samusik et al., 2016) tailored to parsing the multidimen- Given the low cost of converting a scope to this platform sional single-cell data enabled the detection of rare cells types hi hi + (enabled by a simple fluidics device for automated sample (such as CD4 MHCII [Lti] cells, CD11c B cells [ABC cells]) washes in a customized microscope stage), the CODEX tech- and simultaneous placement of those cells in the tissue architec- nique could enable deep studies of various tissue models even ture. Cell-interaction analysis with CODEX recapitulated known with limited resources and instrumentation. features of splenic tissue architecture and revealed that most The unique set of algorithms described here successfully splenic cell types were frequently in homotypic interactions— identified individual cells in the crowded environment of lymphoid which might underscore a novel driving principle of lymphoid tis- tissue by relying on both the information from nuclear and the sue architecture. Further, an important principle was derived from membrane staining. An accurate quantification of single-cell the data-driven i-niche analysis; i.e., we have established that expression data was obtained directly from the images by certain markers, such as B220, CD79b, or CD27, exhibited signif- creating a special algorithm for positional spill compensation. As icant changes in expression levels depending on the tissue of today, this algorithm is only applicable to uniformly distributed context in which the cells reside. As clearly observed in experi- surface markers. Future changes might be required to accommo- ments that drove Figures 4 and 6, cell populations that would date markers that follow a different distribution, i.e., localized to otherwise be thought of as ‘‘broadly’’ expressing a given marker lipid rafts or immune synapses. Nevertheless, the use of this algo- set(Figure 4H), in fact, were composed of multiple subphenotypes rithm enabled us to extract FACS-like data from tissue imaging that correlated with the i-niche identity. In other words, what im- and leveraged the automated phenotype mapping framework munologists previously thought of as a single cell type could be previously developed for CyTOF and multicolor FACS. subdivided into more subtle cell subsets that are defined by the Performance of CODEX on tissue sections was validated in neighborhood in which they reside. We leave open the question analysis of spleen sections of normal and lupus afflicted mice of whether the cells with different properties are attracted to a Cell 174, 968–981, August 9, 2018 979 d EXPERIMENTAL METHOD DETAILS set of neighbor cells, or a given expression level of markers at- B Oligonucleotide sequences tracts the neighbors, or some dialectics thereof. What is clear, B Primer dependent panels to extend the multiplexing however, is that there are more subtle phenotypes in tissues capacity of CODEX than previously assumed, and that future developments of the B CyTOF CODEX comparison technology and the algorithms will shed more light on these B Antibody conjugation, staining and CODEX rendering phenomena. B Primer dependent panels CODEX enabled a quantitative description of autoimmune- B Imaging related changes in the splenic tissue architecture. Among hall- d QUANTIFICATION AND STATISTICAL ANALYSIS marks of MRL/lpr progression were dissipation of marginal B Initial image analysis and segmentation zone, disintegration of PALS, invasion of red pulp with erythro- B Spatial spillover compensation blasts, and the infiltration of mixed-identity B220 DN T cells, B Cell type definition which, interestingly, localize in a niche in between PALS and B Cell interaction analysis the B cell zone and in the marginal zone. A contact-dependent B i-niche analysis effect of B220 DN T cell on CD4 T cells reflected in increased B Neural network training and data analysis levels of activation marker CD27 was observed. An account of d DATA AND SOFTWARE AVAILABILITY statistically significant differences in frequency and strength of pairwise cell-type contacts was created. From these observa- SUPPLEMENTAL INFORMATION tions and their quantitative analysis, we concluded that it is largely the change in cell numbers rather than in cellular interac- SupplementalInformation includesseven figures, fourtables,and two videosand tion strength (estimated from ratio of observed to expected can be found with this article online at https://doi.org/10.1016/j.cell.2018.07.010. probability of interaction) that is involved in reorganization of spleen during transition from norm to autoimmunity. We show ACKNOWLEDGMENTS how i-niche statistics can be used to account for the list of dis- ease-driven changes in sub-splenic anatomy. We also show We thank A. Trejo and A. Jager for technical assistance. We are grateful to Peter Jackson for advice and access to the microscope used for this work. that disease-associated areas of the tissue can be identified This work was further supported by the US NIH, grants and sub awards: independently of the image segmentation, by applying a convo- U19AI057229, 1U19AI100627, R01CA184968, 1R33CA183654-01, lutional neural network to the multidimensional image data, even R33CA183692, 1R01GM10983601, 1R21CA183660, 1R01NS08953301, after training on just one sample. 5UH2AR067676, 1R01CA19665701, R01HL120724, and R01HL128173, Bill Recent advances in genomics suggest that, despite the vast- and Melinda Gates foundation grant OPP1113682, Department of Defense ness of a genetic repertoire, there exist only a limited number of (CDMRP) grant W81XWH-14-1-0180, Northrop-Grumman Corporation agree- cellular states with a concomitantly limited gene expression ment 7500108142, and NIH grant to Stanford Center of Cancer Systems Biology U54-CA209971. G.P.N. is supported by the Rachford & Carlotta A. pattern. These countable, limited patterns are reflected in the Harris Endowed Chair. expression of surface-marker phenotypes recognizable as cell types. It is therefore reasonable to suggest that cell-to-cell AUTHOR CONTRIBUTIONS interactions should be limited as well and falling into repeated patterns. By this token, the data collected in this study lays Conceptualization, Y.G., N.S., and G.P.N.; Methodology, Y.G., N.S., and the foundation for a pan-cellular reference database defining J.K.-D.; Software, N.S. and S. Bhate; Validation, Y.G., N.S., J.K.-D., and G.V.; Formal Analysis, Y.G. and N.S.; Investigation, Y.G. and N.S.; Resources, cellular types not only by identities of proteins expressed, but Y.G., N.S., G.V., S. Black, and M.H.; Data Curation, Y.G. and N.S.; Writing – also by definitions for specific cell-to-cell interactions. We Original Draft, Y.G. and N.S.; Writing – Review and Editing, Y.G., N.S., and performed deep characterization here for normal and diseased G.P.N.; Visualization, Y.G. and N.S.; Supervision, Y.G. and G.P.N.; Project tissue from such a perspective of cell-cell arrangements and Administration, Y.G. and G.P.N.; Funding Acquisition, G.P.N. present here, for the research community, a large (700,000 cells) public dataset encompassing segmentation, quantifica- DECLARATION OF INTERESTS tion, and, most uniquely, spatial data from normal and disease- afflicted spleens (http://welikesharingdata.blob.core.windows. G.P.N., Y.G., and N.S. all have ownership in or are consultants of Akoya Bio, Inc., and members of its scientific advisory board. J.K.-D. was at Stanford net/forshare/index.html). Further analysis of data could enable during the study and currently is an employee of Akoya Bio. Akoya makes re- advances in understanding of clinically relevant cell interactions agents/instruments that are dependent upon licenses from Stanford Univer- in immune tissues as well as development of computational al- sity. Stanford University has been granted a US patent 9909167 that covers gorithms for tissue cytometry and digital pathology. some aspects of the technology described in the paper. Received: September 22, 2017 STAR+METHODS Revised: February 5, 2018 Accepted: July 3, 2018 Detailed methods are provided in the online version of this paper Published: August 2, 2018 and include the following: REFERENCES d KEY RESOURCES TABLE d CONTACT FOR REAGENT AND RESOURCE SHARING Bendall, S.C., Simonds, E.F., Qiu, P., Amir, A.D., Krutzik, P.O., Finck, R., d EXPERIMENTAL MODEL AND SUBJECT DETAILS Bruggner, R.V., Melamed, R., Trejo, A., Ornatsky, O.I., et al. (2011). Single-cell 980 Cell 174, 968–981, August 9, 2018 mass cytometry of differential immune and drug responses across a human he- Parslow, A., Cardona, A., and Bryson-Richardson, R.J. (2014). Sample drift matopoietic continuum. Science 332, 687–696. correction following 4D confocal time-lapse imaging. J. Vis. Exp. Published online April 12, 2014. https://doi.org/10.3791/51086. Burrell, B.E., Warren, K.J., Nakayama, Y., Iwami, D., Brinkman, C.C., and Preibisch, S., Saalfeld, S., and Tomancak, P. (2009). Globally optimal stitching Bromberg, J.S. (2015). Lymph node stromal fiber ER-TR7 modulates CD4+ of tiled 3D microscopic image acquisitions. Bioinformatics 25, 1463–1465. T cell lymph node trafficking and transplant tolerance. Transplantation 99, Robinette, M.L., Fuchs, A., Cortez, V.S., Lee, J.S., Wang, Y., Durum, S.K., Gil- 1119–1125. fillan, S., and Colonna, M.; Immunological Genome Consortium (2015). Tran- Chattopadhyay, P.K., and Roederer, M. (2012). Cytometry: Today’s technol- scriptional programs define molecular characteristics of innate lymphoid cell ogy and tomorrow’s horizons. Methods 57, 251–258. classes and subsets. Nat. Immunol. 16, 306–317. Gabriel, K.R., and Sokal, R.R. (1969). A new statistical approach to geographic Rubtsov, A.V., Rubtsova, K., Fischer, A., Meehan, R.T., Gillis, J.Z., Kappler, variation analysis. Syst. Zool. 18, 259. J.W., and Marrack, P. (2011). Toll-like receptor 7 (TLR7)-driven accumulation of a novel CD11c B-cell population is important for the development of auto- Gerdes, M.J., Sevinsky, C.J., Sood, A., Adak, S., Bello, M.O., Bordwell, A., immunity. Blood 118, 1305–1315. Can, A., Corwin, A., Dinn, S., Filkins, R.J., et al. (2013). Highly multiplexed sin- gle-cell analysis of formalin-fixed, paraffin-embedded cancer tissue. Proc. Rubtsova, K., Rubtsov, A.V., Thurman, J.M., Mennona, J.M., Kappler, J.W., Natl. Acad. Sci. USA 110, 11982–11987. and Marrack, P. (2017). B cells expressing the transcription factor T-bet drive lupus-like autoimmunity. J. Clin. Invest. 127, 1392–1404. Gerner, M.Y., Kastenmuller, W., Ifrim, I., Kabat, J., and Germain, R.N. (2012). Samusik, N., Good, Z., Spitzer, M.H., Davis, K.L., and Nolan, G.P. (2016). Auto- Histo-cytometry: A method for highly multiplex quantitative tissue imaging mated mapping of phenotype space with single-cell data. Nat. Methods 13, analysis applied to dendritic cell subset microanatomy in lymph nodes. Immu- 493–496. nity 37, 364–376. Schubert, W., Bonnekoh, B., Pommer, A.J., Philipsen, L., Bo¨ ckelmann, R., Ma- Jacobson, B.A., Panka, D.J., Nguyen, K.A., Erikson, J., Abbas, A.K., and lykh, Y., Gollnick, H., Friedenberger, M., Bode, M., and Dress, A.W. (2006). Marshak-Rothstein, A. (1995). Anatomy of autoantibody production: Dominant Analyzing proteome topology and function by automated multidimensional localization of antibody-producing cells to T cell zones in Fas-deficient mice. fluorescence microscopy. Nat. Biotechnol. 24, 1270–1278. Immunity 3, 509–519. Socolovsky, M., Murrell, M., Liu, Y., Pop, R., Porpiglia, E., and Levchenko, A. Kanauchi, H., Furukawa, F., and Imamura, S. (1991). Characterization of cuta- (2007). Negative autoregulation by FAS mediates robust fetal erythropoiesis. neous infiltrates in MRL/lpr mice monitored from onset to the full development PLoS Biol. 5, e252. of lupus erythematosus-like skin lesions. J. Invest. Dermatol. 96, 478–483. Spitzer, M.H., Gherardini, P.F., Fragiadakis, G.K., Bhattacharya, N., Yuan, Koh, D.R., Ho, A., Rahemtulla, A., Fung-Leung, W.P., Griesser, H., and Mak, R.T., Hotson, A.N., Finck, R., Carmi, Y., Zunder, E.R., Fantl, W.J., et al. T.W. (1995). Murine lupus in MRL/lpr mice lacking CD4 or CD8 T cells. Eur. (2015). IMMUNOLOGY. An interactive reference framework for modeling a dy- J. Immunol. 25, 2558–2562. namic immune system. Science 349, 1259425. Lieberum, B., and Hartmann, K.U. (1988). Successive changes of the cellular Spitzer, M.H., Carmi, Y., Reticker-Flynn, N.E., Kwek, S.S., Madhireddy, D., composition in lymphoid organs of MRL-Mp/lpr-lpr mice during the develop- Martins, M.M., Gherardini, P.F., Prestwood, T.R., Chabon, J., Bendall, S.C., ment of lymphoproliferative disease as investigated in cryosections. Clin. Im- et al. (2017). Systemic immunity is required for effective cancer immuno- munol. Immunopathol. 46, 421–431. therapy. Cell 168, 487–502. Cell 174, 968–981, August 9, 2018 981 STAR+METHODS KEY RESOURCES TABLE REAGENT or RESOURCE SOURCE IDENTIFIER Antibodies B220 (Ra3-6B2) BD Cat#: 553084; RRID: AB_394614 CD106 (429(MVCAM.A)) BD Cat#: 553330; RRID: AB_394786 CD11b (M1/70) BD Cat#: 553308; RRID: AB_394772 CD11c (HL3) BD Cat#: 117302; RRID: AB_313771 CD16/32 (2.4G2) BD Cat#: 553142; RRID: AB_394657 CD169 (MCA947G) BioRad Cat#: MCA947G; RRID: AB_322322 CD19 (1D3) BD Cat#: 553783; RRID: AB_395047 CD21/35 (7G6) BD Cat#: 553817; RRID: AB_395069 CD27 (LG.3A10) BD Cat#: 124202; RRID: AB_123645 CD3 (17A2) BD Cat#: 555273; RRID: AB_395697 CD31 (MEC13.3) BD Cat#: 553370; RRID: AB_394816 CD35 (8C12) BD Cat#: 558768; RRID: AB_397114 CD4 (RM4-5) BD Cat#: 553043; RRID: AB_394579 CD44 (IM7) BD Cat#: 553131; RRID: AB_394646 CD45 (30-F11) BD Cat#: 553076; RRID: AB_394606 CD5 (53-7.3) BioLegend Cat#: 100602; RRID: AB_312731 CD71 (C2F2) BD Cat#: 553264; RRID: AB_394742 CD79b (HM79b) BioLegend Cat#: 132802; RRID: AB_207588 CD8a (53-6.7) BD Cat#: 553027; RRID: AB_394565 CD90 (G7) BioLegend Cat#: 105202; RRID: AB_313169 ERTR7 (ERTR7) BioRad Cat#: MCA2402; RRID: AB_915429 F4/80 (T45-2342) BD Cat#: 565409 IgD (11-26c.2a) BD Cat#: 553438; RRID: AB_394858 IgM (II/41) BD Cat#: 553405; RRID: AB_394842 Ly6C (HK1.4) Biolegend Cat#: 128001; RRID: AB_113421 Ly6G (1A8) BD Cat#: 551459; RRID: AB_394206 MHCII (M5/114.15.2) BD Cat#: 556999; RRID: AB_396545 Ter119 (TER-119) Biolegend Cat#: 116202; RRID: AB_313703 For composition of oligo conjugated and metal-labeled This paper N/A abs see Table S1 and antibody conjugation part of STAR Methods Biological Samples Spleens from MRL/lpr mice This paper N/A Spleens from BALBc mice This paper N/A Isolated mouse splenocytes This paper N/A Chemicals, Peptides, and Recombinant Proteins dUTP-ss-Cy5; dCTP-ss-Cy3 Jena Custom order Exo- fragment of Klenow DNA polymerase NEB Cat#: M0212L TCEP SIGMA-Aldrich Cat#: C4706-10G Iodoacetamide SIGMA-Aldrich Cat#: I1149-5G Experimental Models: Organisms/Strains MRL/lpr mice Jackson Laboratory IMSR Cat# JAX:000485, RRID:IMSR_JAX:000485 BALBc Jackson Laboratory IMSR Cat# JAX:000651, RRID:IMSR_JAX:000651 (Continued on next page) e1 Cell 174, 968–981.e1–e6, August 9, 2018 Continued REAGENT or RESOURCE SOURCE IDENTIFIER Oligonucleotides (see Table S3) N/A N/A Software and Algorithms CODEX Toolkit This paper https://github.com/nolanlab/CODEX VORTEX environment for X-shift clustering This paper https://github.com/nolanlab/VORTEX Microvolution software for image deconvolution Microvolution http://www.microvolution.com Deposited Data Dataset with data tables can be found at Mendeley https://data.mendeley.com/datasets/zjnpwh8m5b/1 Other Primary image data and additional info can be found at Online repository https://welikesharingdata.blob.core.windows.net/ forshare/index.html CONTACT FOR REAGENT AND RESOURCE SHARING Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Prof. Dr. Garry Nolan (gnolan@stanford.edu). EXPERIMENTAL MODEL AND SUBJECT DETAILS 9 months old female MRL/lpr (IMSR Cat# JAX:000485, RRID:IMSR_JAX:000485) (chosen to represent lupus disease at a pronounced splenomegaly stage) and age/sex matched control BALBc mice (IMSR Cat# JAX:000651, RRID:IMSR_JAX:000651) purchased from Jackson Laboratory were used for the study. All animal studies were done in compliance with ethical regulations and procedures set in the Stanford Administrative Panel on Laboratory Animal Care Protocol 15986. In coherence with the primarily technical purpose of the study no animal cohort randomization or investigator blinding to group allocation was performed. EXPERIMENTAL METHOD DETAILS Oligonucleotide sequences Single base extension during CODEX can be achieved by either a ‘‘missing base’’ approach (Figure 1A) or a ‘‘reversible terminator’’ method (see Video S1, part 2). In the case of the ‘‘missing base’’ approach, which was chosen for the experiments outlined in this paper, the top strand of the double-stranded oligonucleotide is covalently bound to the capture agent (in this case, an antibody) and the bottom strand is annealed through hybridization to the top strand. All antibodies contain the same top strand (5 -ATAG CAGTCCAGCCGAACGGTAGCATCTTGCAGAA-3 ) and different bottom strands. The sequence of the bottom strands contains a common region that hybridizes to the top strand (.TTCTGCAAGATGC 0 0 TACCGTTCGGCTGGAddC-3 ) as well as a 5 variable sequence region that serves as the indexing region. As shown in Figure 1A, the overhanging 5 end of the lower strand of the double-stranded oligonucleotide tag (which forms the overhang) is of the general 0 0 0 0 formula 5 -[C/T] [A/G][5 -C /T -3 ] TTCTGCAAGATGCTACCGTTCGGCTGGAddC-3 The first block a short 5-nt stretch of 5 1-4 1-4 n- random C/T composition designed to increase the polymerase residence on the DNA duplex. The second block is a single nucleotide (either G or A) that allows for incorporation of a labeled dNTPs (dU-ss-Cy5 or dC-ss-Cy3, respectively). The third block is the ‘‘index- ing barcode‘‘ that consists of n random-length homopolymer stretches (1-4 nucleobases each) of alternating ‘‘indexing’’ nucleobases dC and dT that serve as a template for extension of the top oligo with unlabeled nucleotides (dATP and dGTP). Here, n specifies the number or extension cycles after which the fluorescent nucleobase will be incorporated into the duplex. Examples of CODEX index- ing barcodes are CCCTCC for n = 3 and CCCTCCTTTCTT for n = 6. The purpose of having the homopolymer stretches of variable length (e.g., CCCTCCTTTCTT) rather than single base (e.g., CTCTCT) is to increase the polymerase extension specificity and prevent misalignment of upper and lower strands of double-stranded oligonucleotide tags. All oligonucleotide sequences used in this study can be found in Table S1. Primer dependent panels to extend the multiplexing capacity of CODEX CODEX operates using an indexed polymerization step that enables precise incorporation of fluorophores into oligonucleotide-Ab conjugates at predetermined cycles. Although consistent performance of a model antigen (CD45) was observed across 15 cycles of CODEX (Figures S1A–S1F), a gradual accumulation of polymerization errors during each cycle could potentially result in Cell 174, 968–981.e1–e6, August 9, 2018 e2 non-cognate rendering, and thus diminished and/or non-specific signals at later index cycles. In addition, the use of long single- stranded oligonucleotides that would enable indexing beyond 15 rounds might be problematic due to non-specific binding events to tissues under study. For the polymerization event to initiate, a 3 hydroxyl is required. Thus, we reasoned that dedicated primers (each containing a distinct initiating sequence with a 3 hydroxyl) could be used to activate distinct subpanels of antibodies (Figure S2A). This would allow design of antibody panels exceeding 30 markers into subpanels, each with a subpanel-specific activation sequence designed 5 to the indexing region. In this design, the antibody attachment linker is terminated with ddC, such that the extension is only possible after a hybridization of a hydroxyl-containing panel-specific activation primer. The feasibility of such multipanel CODEX design and the robustness of CODEX protocol after many cycles and its independence of staining from the cycle number were tested in a model experiment. A 22-color panel of antibodies (11 cycles) conjugated to a termi- st nd rd nated top oligonucleotide, was hybridized with lower oligonucleotides of 1 ,2 , and 3 panels (Figure S2B). Thus, every antigen is detected thrice by the same antibody conjugated to oligonucleotides of 3 different panels. Each panel can only be rendered after annealing of a panel-specific activator oligonucleotide. The staining was rendered in 36 cycles (11 detection cycles + 1 blank no-anti- body cycle per activator oligo) of CODEX with additional activator oligonucleotide hybridization step between each of the 3 panels. st th th The signal for same antibody detected at different cycles (e.g., 1 ,13 , and 24 ) was consistent across the three panels (Figure S2C). This panel-activator design extends CODEX to a theoretically unlimited multiplexing capacity, bounded only by the speed and res- olution of the imaging process itself and the time required for each imaging cycle. CyTOF CODEX comparison Cell preparation and staining by metal tagged antibody for CyTOF analysis was performed as described before (Spitzer et al., 2017). TM Mass cytometry was performed on a CyTOF 2 mass cytometer (Fluidigm) equilibrated with ddH2O. For CODEX analysis, isolated spleen cells were stained a panel of antibodies conjugated to indexing oligonucleotides. Samples were fixed to a coverslip (Figure 2A) and imaged over 12 cycles of CODEX protocol. Images were segmented using the in situ cytometry software toolkit developed for this study (see Figure 2A for exemplary segmentation of the cell spread), and the staining of individual cells across the indexing cycles was quantified. Segmentation data was converted into flow cytometry standard (FCS) format and analyzed using the conventional flow cytometry analysis software Cytobank. Antibody conjugation, staining and CODEX rendering Detailed stepwise CODEX protocols can be found online (see Key Resources Table and Data and Software Availability section below). For full list of antibody clones and vendors see Table S1. Custom manufactured microfluidic setup (Figures S7A–S7C) was used to automate CODEX solution exchange and image acquisition. Instrument and blueprints and control software are available upon request. Primer dependent panels Rendering of antibodies with spacers followed the same procedure as the standard CODEX protocol with the exception of the following differences. Before proceeding to rendering next spacer dependent panel, the stained cells were incubated with a spacer oligonucleotide (1 mM final concentration in buffer 405) at room temperature for 10 minutes. Cells were washed 4X with buffer 4 and rendering proceeded as usual. To initiate each additional spacer set, the spacer incubation step was repeated using corresponding spacer samples. Imaging Images were collected using a Keyence BZ-X710 fluorescent microscope configured with 3 fluorescent channels (FITC, Cy3, Cy5) and equipped with Nikon PlanFluor 40x NA 1.3 oil immersion lens. Imaging and washes were iteratively performed automatically us- ing a specially developed fluidics setup (Figures S7A–S7C). Images were subject to deconvolution using Microvolution software (http://www.microvolution.com/). The staining patterns of 28 DNA-conjugated antibodies were acquired over 14 cycles of CODEX imaging and overlaid with 2 additional fluorescent antibodies, CD45-FITC and NKp46-Pacific Blue and a DRAQ5 nuclear stain (Figure 3A and low-resolution views in Video S2). Each tissue was imaged with a 40x oil immersion objective in a 7x9 tiled acquisition at 1386x1008 pixels per tile and 188 nm/pixel resolution and 11 z-planes per tile (axial resolution 900 nm). Images were subjected to deconvolution to remove out-of-focus light. After drift-compensation and stitching, we obtained a total of 9 images (one per tissue) with x = 9702 y = 9072 z = 11 dimensions, each consisting of 31 channels (30 antibodies and 1 nuclear stain). QUANTIFICATION AND STATISTICAL ANALYSIS Initial image analysis and segmentation For each imaging field analyzed by CODEX multidimensional staining multi-color z stacks collected during individual cycles were aligned against reference channel (CD45) by 3D drift compensation (Parslow et al., 2014). If necessary individual fields covering large tiled areas were ‘‘stitched’’ using dedicated ImageJ plugin (Preibisch et al., 2009). For the 22-channel experiment on dissociated cells attached to coverslip (Figure 2) images corresponding to the best focal plane of vertical image stacks collected at each acquisition e3 Cell 174, 968–981.e1–e6, August 9, 2018 step of CODEX were chosen for quantification. For the 31-channel main experiments on mouse spleen sections, the segmentation was performed on the whole image stack using a volumetric segmentation algorithm described below. For this study, we purposefully developed a 3D image segmentation that combines information from the nuclear staining and a ubiquitous membrane marker (in this case CD45) to define single-cell boundaries in crowded images such as lymphoid tissues. This algorithm inverts the membrane image and multiplies it with the nuclear image, creating a synthetic image with enhanced contrast between neighboring nuclei. This image is subject to low-pass FFT filtering and an then individual cell objects (collections of voxels) are identified using a gradient-tracing water- shed algorithm. Per-cell intensities were quantified by integrating the intensity of each channel within a given cell object and divided by the region size in voxels. We benchmarked the segmentation algorithm against a dataset of BALBc mouse spleen images with expert hand-labeled nuclei in and we found the algorithm was able to correctly identify 87.25% ± 2.89% cells, of which 89.88% ± 1.12% were singlets (one-to-one correspondence between a hand-labeled cell and a segmented object) (Figures S1H–S1J). For each segmented object (i.e., cell) a marker expression profile, as well as the identities of the nearby neighbors were recorded using Delaunay triangulation (https://data. mendeley.com/datasets/zjnpwh8m5b/1). Spatial spillover compensation Accurate segmentation per se is not sufficient to obtain high-quality estimate of single-cell expression from images. The reason for that is that in lymphoid tissue the cells are so tightly adjacent that their membrane signals can partially overlap, resulting in blending of signals between neighboring cells, the phenomenon termed spatial spillover. In order to compensate for that, we estimated the cell- to-cell signal spill coefficients based on the fraction of shared boundary between each pair of cell objects, resulting in a banded matrix (most cells don’t have any shared boundaries). To compensate the cell-to-cell spill, the raw intensity vector was multiplied by the inverse spill matrix (Figure 2D). Besides the spatial signal spillover, there are other factors that add to artifactual cell-like objects: debris misidentified as cells, dou- blets (two adjacent cells merged together) as well as autofluorescent objects, both of which can lead to spurious as double-positive signals on the biaxial scatterplots. By analogy to how debris and doublets are eliminated from FACS data by applying special ‘Singlet’ gates to SSC-FSC parameters, we devised a ‘cleanup’ gating strategy based on several quality control parameters: nuclear stain density (nuclear signal divided by cell size), profile homogeneity (relative variance of signal from cycle to cycle), background staining on blank cycles and, finally, nuclear signal and cell size (Figure S1K). We found that applying those filtering gates had a synergetic effect with the compensation, reducing the frequency of spurious double-positive cell signals by approximately an additional factor of 2 (Figure S1L, e.g., compare fraction of CD5/IgD double positive cells in Ungated-Compensated and Post Cleanup gated – Compensated in (L)). Cell type definition The 9-spleen dataset was subject segmentation, quantification, compensation and cleanup gating, as described above, yielding a total of 734101 30-dimensional single-cell protein marker expression profiles (Figure 3C, https://data.mendeley.com/datasets/ zjnpwh8m5b/1). The segmented CODEX data was subject to automated phenotype mapping algorithm X-shift that was previously developed and validated on CyTOF data (Samusik et al., 2016)(Figure 3C). 58 phenotypic clusters inferred by X-shift clustering were manually annotated (Figures 3C and 3D and Table S2) based on the 30-color marker expression profile and thorough visual inspec- tion of the representative image samples (Figures S5A–S5J, more examples of ‘‘cell passports’’ can be found in associated online repository - see STAR Methods). Some clusters were found to originate from imaging artifacts such as dust and tissue sectioning defects. That reduced the overall number of cell-like objects to 707466. Each cluster was assigned to one of 27 broadly defined sin- gle-cell phenotypic groups (cell types), which in some cases could be clearly matched to major immune cell types and in others were named according to expression of distinguishing surface markers (see cluster annotation and cell counts in Table S2). Cell interaction analysis To define for each cell the neighbors of the first (immediate) tier of proximity Delaunay graph was computed for the dataset (https:// data.mendeley.com/datasets/zjnpwh8m5b/1). The odds ratio of co-occurrence of cell type A and cell type B was estimated as the observed frequency of co-occurrence (mean of the beta-distribution, with parameter alpha = number of edges connecting cell types A and B and parameter beta = total number of edges minus number of edges connecting A-B) divided by the theoretical frequency of co-occurrence (total frequency of edges incident to type A multiplied by the total frequency of edges incident to type B) see Table S3. The odds ratios are represented in heatmaps on Figure 3G, with a range of values from less than 1 to more than 1 mean- ing that two cell types are, respectively, less or more likely to co-occur than expected by chance. The significance of the difference from zero was tested using binomial distribution (probability of getting an observed number of interactions between A and B (successes) among the total number of registered interactions (number of trials) given the theoretical probability of A-B interaction (probability of success)). The significance of change of interaction frequencies or log-odds ratios were computed between BALB/c and Stage 1 (early) MRL using pairwise t test. However, the same procedure could not be applied to testing BALB/c versus MRL/lpr Stages 2 and 3 because of high sample-specific variation in those more advanced disease stages. Therefore we scored computed the deviation of those Stage2/3 values from BALB/c using c statistics because it does not require Stage 2/3 samples to have a common mean. Cell 174, 968–981.e1–e6, August 9, 2018 e4 The P values were subject to FDR correction using Benjamini–Hochberg procedure. Interactions that were considered significant for FDR q-value < 0.05 or > 0.95 (Table S3). In order to estimate the overall deviation of either interaction frequency matrices or log-odds ratio matrices, the matrices were sub- ject to z-transformation based on the mean and the SDs of the BALB/c samples, and then c statistics was computed as square root of the sum of squares of all elements of the z-score transformed matrices (Figure 5F). i-niche analysis The i-niche analysis was performed based on 2-dimensional Delaunay triangulation of the cell center coordinates. Delaunay triangu- lation and the related concepts of Voronoi Tesselation and Gabriel graphs were previously applied in eco-geographical analysis of species distribution (Gabriel and Sokal, 1969) and therefore were deemed as equally applicable to the analysis of tissue organization on the single-cell level. The i-niche is defined as a set of first-order Delaunay neighbors of the given ‘index cell’, i.e., the i-niche cells are the ones that are directly connected to the i-cell with edges in the Delaunay triangulation of cell centers. We distinguish i-niche from the more formal understanding of ‘‘niche,’’ which is often used in stem cell literature and where numbers of cells in the niche and their placement within the niche is undefined. In our definition, we allow the central cell to be of any type and are counting the cell types present in the ring. This flexible definition allows for multi-cellular interactions around a central cell to define the biology of that cell (and vice versa). Computationally, the i-niche window slides from cell to cell, considering each set of adjoining cells— and therefore allows consideration of the constituencies of different central cell types that might populate a given i-niche. We under- stand that our current definition is arbitrary and could be extended to include other specific cell arrangements—including, though beyond the scope of the current work, a 3D sphere of cells contacting the index cell. Neural network training and data analysis Preprocessing Image stacks were maximum-intensity projected following deconvolution. Data was quantile normalized to 4 levels (0, 0.25, 0.5 and 0.75 quantiles). A baseline model was able to distinguish models without this discretization and normalization, suggesting strain-spe- cific differences in antibody staining intensity. Training and cross validation split Four spleen samples (two BALB/c and two MRL/lpr) were chosen as training samples. The remaining five spleens tissue samples (one BALB/c and four MRL/lpr) were used for testing the trained model. For cross-validation, different combinations of spleens were allo- cated to training and test sets. During training, 224x224 images were randomly extracted from the training tissue samples, at 1x, 0.5x and 2x zoom. At 1x zoom, there would be 6804 non-overlapping image patches in the training dataset. The trained models were tested on 4500 patches, at 1x zoom. Hyperparameters were manually tuned on 500 randomly selected images from the testing spleens. The Adam optimizer was used for training with an initial learning rate of 0.0001. Baseline model. A logistic regression model was trained by averaging marker intensities across the image. L2 regularization was used for weights. Neural network architecture. To avoid the learning of trivial sample-specific staining variation, data were quantile normalized sam- ple-wise and each marker was discretized to four levels. Since disease-specific hallmarks could potentially be present at multiple scales, the training data for our neural network was extracted at multiple levels of magnification. A simple regularized logistic regres- sion model that considered only average marker expression and did not incorporate spatial information was unable to successfully distinguish patches normal and MRL/lpr spleens, whereas the trained neural network model consistently achieved a 90% precision of classification of image patches during cross-validation. A fully convolutional network architecture was used, with the following layers. To generate a prediction for an entire image patch, a global max-pooling layer was used. 1. Conv3 60 2. Conv3 120 3. Conv3 64 4. Batch Norm 5. Conv3, 64 6. Max pooling 2x2 7. Conv3, 128 8. Conv3, 128 9. Max pooling 2x2 10. Conv3, 256, 11. Conv3,256 12. Conv3,256 13. Max pooling 2x2 14. Conv3,512 15. Conv3,512 e5 Cell 174, 968–981.e1–e6, August 9, 2018 16. Conv3,512 17. Conv1,256 18. Conv1,64 19. Conv1,1 20. Global max pooling 21. Sigmoid Weights for layers 5-16 were initialized from the VGG-16 pretrained model. The model was trained with cross-entropy loss. Regularization. L2 regularization (0.1) was used for network weights. L1 regularization was applied to the feature map output after layer 19 to encourage sparse activations Whole sample activations for test set. Since the network was fully convolutional, it could be applied to images of any dimension. The network was applied to entire fields of view individually. The activation maps were obtained as the output after layer 21. Aligning cell type information. Each cell was assigned the MRL/lpr score of the corresponding pixel in the image. Enrichment and neighborhood analysis. FDR controlled chi-square tests of proportions were carried out to determine enrichment of specific cell types in the top 10% of cells by MRL/lpr score. For neighborhood analysis of dendritic cells, the composition of the neighborhoods (cell centers within 30 pixels) of the top 300 cells (by MRL/lpr score) were compared to the composition of the neighborhoods of the bottom 300 cells. Only cells with positive neural network assigned MRL/lpr score, in MRL/lpr regions, were considered for this analysis. DATA AND SOFTWARE AVAILABILITY Software used in the paper for parsing image data can be obtained at: https://github.com/nolanlab/CODEX Data tables can be downloaded from Mendeley: https://data.mendeley.com/datasets/zjnpwh8m5b/1 All primary image data, high resolution focused montages, complete single cell data tables and various additional information can be obtained at: http://welikesharingdata.blob.core.windows.net/forshare/index.html Flow formatted segmented data can obtained from online repository page (link above) or from Cytobank: CODEX on spreads of isolated mouse splenocytes (Figure 2B): https://community.cytobank.org/cytobank/experiments/69534 https://flowrepository.org/experiments/1686 CyTOF on isolated splenocytes (Figure 2B): https://community.cytobank.org/cytobank/experiments/69533 https://flowrepository.org/experiments/1687 CODEX on BALBc spleen tissue sections (Figure 2E): https://community.cytobank.org/cytobank/experiments/69889 https://flowrepository.org/experiments/1688 Cell 174, 968–981.e1–e6, August 9, 2018 e6 Supplemental Figures Figure S1. Benchmarking CODEX, Related to Figures 1 and 2 (A) Experimental scheme for mimicking the tissue with 30 distinct cell types. (B) Montage of a fragment of imaging field of the 15 cycles of CODEX used to render the mix of 30 barcoded spleens – first cycle top left last cycle bottom right. (C) Heatmap (cycles in columns, cells in rows) showing mean fluorescence per cell membrane for each cell per in each of the 15 CODEX cycles performed on cells of 30 barcoded spleens. Odd columns correspond to imaging after labeled base incorporation. Even columns correspond to imaging after inactivation of staining by TCEP. (D) Time-lapse profile of median intensity per cell membrane for individual cells marked by white arrows on (B). (E) Average intensity of CD45 antigen expression in ‘‘positive’’ (blue columns) and ‘‘negative’’ (red columns) cells in 15 CODEX cycles of the experiment. (Similar results were obtained for Cy3-positive populations – data not shown). Linear regression was performed to indicate trends in accumulation of background and signal decline associated with cycle number. (F) Table summarizing CODEX performance stats. Average signal to noise ratio was estimated from ratio of average signal of all positive cells across all cycles to the signal of all ‘‘negative’’ cells across all cycles. Efficiency of fluorophore removal was estimated from average ratio of ([signal after TCEP in cycle N]- [signal after TCEP in cycle (N-1)])/[signal in cycle N] for cells positive in cycle N across all cycles. Average expected signal deterioration was estimated using the trendline equation from (E). Average background accumulation was estimated by fitting linear trendline into the per cycle ratio of average background to average signal (not shown). (G) Image quantification approach used on CODEX data from (A): best focal planes of CODEX stacks were segmented by Cell Profiler. To account for local background the value corresponding to difference between the mean intensity value inside ‘‘cell membrane’’ object (left panel) and the mean intensity inside the (legend continued on next page) external ring object (right panel) was chosen as a representation of the intensity of the antibody signal. In all other experiments custom (see STAR Methods) segmantation developed in this study was used. (H) Sample 500x500 px regions from two samples (BALBc-3 and MRL-4) showing hand-labeled cell centers (yellow crosshairs) and cell outlines detected by the segmentation algorithm (randomly colored). (I) Comparison between the hand-labeled cell identification and algorithm-based algorithm identification, expressed in 3 measures: %Nuclei found (how many of the hand-labeled nuclei centers ended up inside the segmented regions), % Singlets (how many of the cell regions with at least one hand-labeled nuclei center contained exactly one cell center) and %Unlabelled regions (how many segmented regions did not contain a hand-labeled cell center). (J) Summary statistics comparing the segmentation quality between BALBc and MRL/lpr samples. (K) Three step cleanup gating strategy based on 1) stain density (nuclear signal divided by cell size) and profile homogeneity (relative variance of signal from cycle to cycle), 2) removing objects with high background by gating on the signal accrued in ‘‘blank’’(no stain) cycles 3) constraining the cell size. (L) Percentage of artifactual double-positive cells in CODEX data from sample BALBc-2 (as seen in the upper right quadrant biaxial flow style plots of mutually exclusive lineage markers IgD and CD5) depending on gating and spill compensation. Figure S2. Expanding the Multiplexing Limit of CODEX by ‘‘Panels and Activators’’ Design, Related to Figure 1A and STAR Methods (A) Diagram of ‘‘multipanel’’/’’activator oligo’’ CODEX approach. The list of antibodies can be divided in sets such that number of antibodies in each individual set does not exceed the capacity of the multiplexing protocol to render staining without significant signal loss (e.g.30). Each such set of antibodies will be conjugated to ‘‘terminated’’ (the last 3 base is dideoxy- or propyl- modified) upper strand oligonucleotide of the same sequence as in the original version of the ‘‘missing base’’ approach. The lower strand oligonucleotides will incorporate an additional set-specific region, which will serve as a landing spot for the dedicated primer oligo which is to be on-slide hybridized to the particular subset of the total plurality of the antibodies at the time when they are to be rendered. This approach prevents extension of reads beyond certain threshold and at the same time have an unlimited potential number of antibodies in the sample. (B) Schematics of experiment demonstrating the ‘‘activator’’ method and its robustness. Each antigen of a set of 22 surface markers is redundantly detected by three CODEX tag conjugates of the same antibody. The first conjugate is detected during panel 1 rendering, second – during panel 2 etc.. Thus the signal for same st th th antigen is detected at different cycles (e.g., 1 ,13 , and 24 ). (C) Montage of a fragment of imaging field of the 36 cycles of CODEX used to render a mixture of 18 barcoded spleens (similar to design in Figure 2A). Cycles N,N+12 and N+24 all three of which render same pair of antigens are shown per tile for all 11 pairs of antigens (see annotation in the black rectangle of each tile). Figure S3. Types of Samples in MRL/lpr Dataset, Related to Figure 5 MRL/lpr dataset has 9 samples: 3 control wild-type BALBc spleens (BALBc 1,-2,-3 and 6 MRL spleens MRL 4,-5,-6,-7,-8,-9). Based on disintegration of marginal zone as measured by frequency of marginal zone macrophages (MZM’s, – see black asterisk on Figure 5B and yellow arrow in this figure pointing to the area where CD169 positive (red) rim of MZMs is expected to be observed) and accumulation of double negative T cells expressing B220 B cell marker (B220 DN T cells – see red asterisk on Figure 5B) MRL spleens were grouped into early (MRL 4,-5,-6), intermediate (MRL 7,-8), and late (MRL 9) types. Early stage was represented by 3 MZM positive DN T cell-low spleens. Two spleens represented the intermediate stage: MZM low DN T cell-low spleen (Int1) and MZM positive DN T cell-positive spleen (Int2). Late stage was represented by single MZM positive DN T cell-positive spleen. A single representative spleen is shown for each stage together with interaction matrix. Color represents odd ratios (observed frequency of interaction/ expected frequency of interaction). Figure S4. CODEX Pinpoints Splenic Location of Unique Cell Types, Related to Figures 3–5 (A) Distribution of CD4(+)MHCII(+) cells (marked with white circles) in BALBc #2 spleen stained with IgD (green) and CD90 (red) to indicate positions of B and T cells accordingly. (B) CD4 and MHCII expression in isolated mouse splenocytes gated negative for all CODEX panel markers and in addition 120 g8 (lineage depletion with BD 558451 and dump channel for FITC conjugated or biotinilated antibodies corresponding to the antigens stained with CODEX panel were used for negative gating) except CD4, MHCII, CD45 and CD44. (C) CD4(+)MHCII(+) cells within the gate shown in (B) were sorted out and subjected to microarray analysis. CD4 T cells, CD8 T cells, bulk B cells and Conventional CD11c positive dendritic cells were co–sorted as a control. Expression of Lti signature genes (two individual signature sets as inferred in (Robinette et al., 2015)) in sorted cells. (D and E) CD11c+ B cells (age associated B cells (ABCs) in normal nd M/lpr spleens. ABCs have been shown to be a key participant in the triggering of certain autoimmune responses (Rubtsova et al., 2017, Rubtsov et al., 2011)) their splenic location has not been previously described in the literature. We observed ABCs to tightly associate with conventional dendritic cells (cDC) and occupy a distinct peri-follicular space in the boundary between PALS and B-zone. Interestingly, these cells diminished in numbers and redistributed toward intra-follicular space in the MRL/lpr spleens. (legend continued on next page) (F and G) Co-distribution of B220 and TCRb in isolated splenocytes of normal (BALBc) and autoimmune (MRL/lpr) mice. Gate in (G) points to significant (13%) presence of B220+ DN T cells in MRL spleen. (H–J) Thread like arrangement of CD8 T cells (purple, annotated with V-letter) has been noticed in PALS of splenic samples across dataset. To examine potential mechanisms driving these structures CD8 Tells and B220 positive B cells were sorted individually from BALBc spleen (I) and later combined in flat bottom microwell plates and mixed at 37C in culture medium. After mixing cells were stained for B220 (green) and CD8a (red) and imaged (J). Thread like structures similar to what was observed in spleen were detected. (K) Heatmap showing average frequencies of cell types (rows of heatmap) in the ring of index cell neighbors (see schematics on the right) for all niche clusters (0-99 in columns). (L) Heatmap shows how different cell types (in rows) are distributed between niches (in columns). Figure S5. ‘‘Cell Passports’’ of Selected Cell Types Identified in Normal and MRL Spleens, Related to Figures 3F and 5B (A) Diagram of per cycle markers for CODEX cycle montages in B,C and D. (B, E, and H) High resolution montage of CODEX cycles with cells of interest (CD11c(+) B cells) marked with yellow crosses is shown in (B). Low resolution montage of distribution of cells of interest (marked with white circles) in all imaged samples is shown in (E). Average expression profile of all markers in the cells of the selected cell type is shown in (H). (C, F, and I) Same for CD4(+)MHCII(+) cells. (D, G, and J) Same for CD106(+)CD16/32(-)Ly6C(+)CD31(+) cells. More examples of ‘‘cell passports’’ can be found in associated online repository (see STAR Methods). Figure S6. Cross-Tissue and Cross-Samples Distribution of Interacting Cell Pairs for Selected Types of Cell-to-Cell Interactions, Related to Figure 5C Interacting cell pairs are marked with white and cyan circles on the montage of IgD CD90 (B cell and T cell markers) staining of every sample of the dataset. Due to cell proximity in most cases cyan circles practically completely overlay white. (A and B) Interaction of CD4 and CD8 T cells with ERTR stroma (change in odds ratio score correlates with change in interaction count). (C–E) Interaction of granulocytes with CD4 T cells, dendritic cells and erythroblasts. (F and G) Interaction of erythroblasts with stromal and B220(+) DN T cells. Interactions in (C-G) scored as increased in early MRL/lpr (4,-5,-6) as compared to BALBc spleens (FDR of t test on normalized interaction counts between conditions < 0.05, difference in interaction counts > 0). (H and I) Interactions of B220(+) DN T cells with CD8 T cells and stromal cells. These interactions scored as increased in intermediate and late MRL/lpr (7,-8,-9) as compared to early MRL/lpr spleens (FDR of t test on normalized interaction counts between conditions < 0.05, difference in interaction counts > 0). More examples of cell type pairs with change in interactions across dataset can be found in associated online repository (see STAR Methods). Figure S7. The Fluidic Setup and Stage for Running CODEX Experiments, Related to STAR Methods (A) General diagram of robotic fluidic setup used in this study. CODEX experiments are done in an open flow cell, which can be imaged in any inverted mi- croscope. Six solutions have to be programmatically delivered and removed from the flow cell, which in the meantime sits in spatially defined position in the imaging system. A combination of 6-channel Tecan syringe pump equipped with 250ul syringes and USB-relay driven vacuum valve was used for iterative solution delivery and removal. Imaging was performed in Keyence BZ-X710 fluorescent microscope configured with 3 fluorescent channels (FITC, Cy3, Cy5) and equipped with Nikon PlanFluor 40x NA 1.3 oil immersion lens. Insets show photographs of actual microscope stage and fluidics robot. (B) Detailed 3D model of CODEX stage used in experiments. A metal insert was machined to be compatible with either ASI (Advanced Scientific Instrumentation) or Keyence 3d stages. Disposable (one per experiment) acrylic platform with a circular cutout in the middle was custom designed and lasercut such that it could be attached to the metal stage insert. Before multicycle run the coverslip with a sample was glued to the acrylic base which produced an open flow cell. As opposed to closed, open flow cell design ensures efficient (99.9%) and rapid solution exchange that is critical for CODEX protocol. (C) An exemplary photograph of full CODEX setup when attached to an inverted confocal microscope.

Journal

CellUnpaywall

Published: Aug 1, 2018

There are no references for this article.