Monkeybread enables investigations into how cells organize into cellular niches, here defined as regions within the tissue microenvironment comprising a specific composition of cell types or phenotypes. To do so, Monkeybread first selects the neighboring cells of each cell (those that fall within a pre-specified radius of each cell) and counts the number of cells of each cell type or phenotype to form an ×

neighborhoods matrix, , where is the number of cells and is the number of cell type or phenotype labels. Element , represents the number of neighbors of cell that are of cell type . To account for neighborhoods with differing numbers of cells, Monkeybread provides the option of normalizing these neighborhood counts by dividing each cell9s neighbor counts by its total number of neighbors. To account for the large differences in cell abundances between cell types, Monkeybread also provides the option of computing the z-scores of each neighborhood count (or normalized count), thus placing count values between cell types on similar scales. Finally, niches are computed by clustering the rows (corresponding to cells) of the normalized neighborhoods matrix. By default, Monkeybread uses the Leiden clustering algorithm 22.

Monkeybread provides a lightweight statistical test for testing enrichment of cellular colocalization between pairs of cell types (Figure S1). This test considers the following test statistic: for two cell types and under test, for each cell of type , we find the distance to its nearest cell of type . We then count the number of such shortest distances from cells of type to type that fall below some user-provided threshold distance . To test for significance, a permutation test is performed that generates a null distribution of the test statistic by permuting the spatial coordinates of nontype cells. This test can optionally be run on cells within specific niches to condition the analysis on co-localization patterns explained by the organization of the tissue into cellular niches. When visualizing the full tissue specimen, there often exist so many cells that it becomes difficult to visually discern the spatial distribution of cells of a given cell type, especially rare cell types, or those expressing a given gene. Monkeybread provides functionality for visualizing and comparing the spatial density of cells with a specific cell type label or phenotype across the tissue. More specifically, Monkeybread uses kernel density estimation to visualize regions with a high concentration of cells exhibiting rare or specific phenotypes (Methods). In addition, Monkeybread provides functions built atop scanpy for easily zooming in on specific regions of the tissue and sub-setting the data to the zoomed-in regions within the smaller field of view.

Lastly, Monkeybread implements a lightweight statistical test developed by He et al. (2021) 23 for identifying intercellular communication via ligand-receptor binding (Methods; Figure S2). This method is related to a spatial ligand-receptor co-expression test implemented in Giotto but here is implemented in Python. It differs from that implemented by Squidpy in that Squidpy9s uses CellPhoneDB 24, which does not consider spatial proximity between cell types. In addition to this test, Monkeybread provides a visualization tool, built atop matplotlib 25 and scanpy 19, that depicts which pairs of cells in the tissue may be communicating via a given ligand-receptor by drawing a line between each pair of neighboring cells of types and that are expressing the given ligand and receptor. These lines are colored according to the magnitude of their ligand-receptor coexpression scores. This visualization tool enables the user to quickly see where in the tissue neighboring cells are co-expressing the ligand-receptor pair.

We applied Monkeybread to a human melanoma sample from Vizgen9s immuno-oncology public data release 21 to study the cellular organization of the TME. First, we annotated the cells with both a coarse and granular set of cell type labels by first clustering the cells and manually crossreferencing the genes expressed in each cluster with known marker genes (Figure 1a,b, S3-7; Methods). We examined the spatial distribution of the broad cell types using Monkeybread9s tool for visualizing cell spatial densities and found malignant cells occupying the core of the tumor surrounded by regions of dense stromal and immune cells (Figure 1.c).

Within the monocyte/macrophage compartment, we identified three distinct populations (Figure 1d). The first displayed a putatively immunogenic phenotype as indicated by high expression of CXCL9, CXCL10, and CXCL11 26. The second displayed markers commonly associated with a tolerogenic, M2-like phenotype (CD163 and MRC1) as well as other genes associated with immune suppression (C1QC and IL10) 27,28. A final population expressed SPP1, which has also been associated with immune suppression and poor clinical outcome (Figure 1.e) 29.

To further understand the immunological role of these macrophage populations, we compared their transcriptional profiles to macrophages differentiated from PBMC-derived monocytes in vitro under three conditions: M1-like macrophages were developed in vitro through polarization of human monocytes by IFNg and TLR4 stimulation. This population is considered to have anti-tumor functions including release of proinflammatory cytokines, phagocytosis, and ability to crosspresent extracellular antigens produced by dying tumor cells to activate CD8 T cells 30. Second, M2-like macrophages were developed in vitro through intensive IL4 stimulation of monocytes. M2 macrophages are known to promote many pro-tumorigenic processes such as angiogenesis, immune suppression, hypoxia induction, tumor cell proliferation, and metastasis 31. Lastly, we deployed a third population that represents a better in vitro model for tumor associated macrophages (TAM), developed under conditions more closely resembling the TME 32. We performed bulk RNA-seq on these three populations and performed differential expression (DE) analysis between the M1-derived population and the M2-derived population as well as between the M1-derived and TAM-derived populations (Methods; Table S1). We found 134 genes higher in M1 vs. M2 and 144 genes higher in M1 vs. TAM (aj.dp < 0.05). We took the intersection of these two gene sets to derive an M1-signature consisting of 101 genes that we used to score the in vivo macrophages (Methods). Similarly, we found 80 genes higher in the M2 vs. M1 populations and 78 genes higher in the TAM vs. M1 populations (adj. p < 0.05). We took the union of these two genes sets to derive an M2/TAM-signature comprising 100 genes that we used to score the in vivo macrophages. Overall, we observed high expression of the M1-signature in the CXCL9high macrophages and high M2 and TAM signatures in the CD163-high and SPP1-high macrophages (Figure 1.f,g). In summary, these results support the immunogenic function of the CXCL9-high and tolerogenic function of the CD163-high and SPP1-high macrophages in the tumor microenvironment.

Next, we examined the spatial distribution of these three macrophage populations using Monkeybread9s kernel density estimation tool and found that these cells occupy distinct regions within the tissue (Figure 1.h). Interestingly, the CXCL9-high and CD163-high populations were concentrated in regions with higher immune cell abundance (immune <hot= regions) and in regions surrounding the tumor. In contrast, the SPP1-high macrophages were found in regions with fewer immune cells (immune <cold= regions) and mostly found within the tumor in regions of dense malignant cells. Although the gene expression signatures of the CD163-high and SPP1-high macrophages indicate a suppressive function, their distinct spatial distribution indicates that they are playing further differentiated roles in the TME.

We then used Monkeybread to identify distinct cellular niches of immune cells in the tumor microenvironment (Figure 2.a). These niches included both immunologically cold (Niches 1 and 10) and hot niches (all remaining niches; Figure 2.b). Niche 7 was enriched for CXCL9-high macrophages, Niches 0 and 10 for CD163-high macrophages, and Niche 1 for SPP1-high macrophages. Niche 7 contained a higher enrichment of diverse immune cells than Niches 0, 1, and 10. We then examined the mean expression in immune cells within each niche of a ten-gene interferon-gamma related signature that has been shown to correlate with response to immunotherapy 33 (Fig. 2.c). Niche 7 displayed high overall expression of genes in this signature while Niche 10 displayed markedly low expression. Moreover, Niches 0 and 10 (the CD163-high enriched niches) displayed marked lower levels of IFNG than the other niches. Together, these results indicate that Niche 7, and in turn CXCL9-high macrophages, have an immunogenic role whereas CD163-high macrophages may have a suppressive role in the TME.

Next, we used Monkeybread to plot the distribution of the number of T cells that are neighbors of CXCL9-high, CD163-high, and SPP1-high macrophages. Consistent with the niche analysis, we found that the CXCL9-high macrophages are surrounded by the most T cells whereas SPP1-high macrophages are surrounded by the fewest (Figure 2.d). Next, we ran Monkeybread9s cell colocalization statistical test between all pairs of myeloid cell subtypes and all T cell subtypes. We found that CD163-high and SPP1-high macrophages only displayed enriched co-localization with T regulatory (Treg) cells and PDCD1-high CD8 T cells whereas CXCL9-high macrophages also colocalized with T follicular helper cells (Tfh) cells and unclassified CD8 T cells (Figure 2.e). These results were consistent with the niche analysis, which showed higher Tfh enrichment in Niche 7 than in Niches 0, 1, or 10. This is a highly relevant observation since Tfh cells have been shown to correlate with reduced tumor growth and improved survival in multiple solid tumors types 34. Together these results further support that these macrophage populations play differing immunological roles within distinct niches in the TME.

Lastly, we used Monkeybread to test for co-expression of the ligand receptor pairs CXCL9/CXCR3 and CXCL10/CXCR3 between all macrophages and their neighboring conventional T cells and found significant co-expression of CXCL9 and CXCR3 (p = 0.001) as well as between CXCL10 and CXCR3 (p=0.040). Using Monkeybread9s function to visualize this co-expression, we confirmed higher co-expression of CXCL9/CXCR3 in Niche 7, which is enriched for CXCL9-high macrophages (Figure 2.f,g), compared to Niche 0, which is enriched for CD163-high macrophages (Figure 2.h,i).

All authors were paid employees of Immunitas Therapeutics, Inc. at the time this work was Monkeybread uses kernel density estimation to visualize the spatial distribution of cells that exhibit a certain phenotype of interest (e.g., T cells), which smooths out the discrete cell locations enabling easier visualization of regions with high or low density 39. Specifically, let * =! denote the spatial coordinates of a cell exhibiting the target phenotype. We assume that there exists an underlying probability distribution ( ) over that we seek to estimate via kernel density estimation. To visualize the concentration of cells across the tissue, we seek to color each cell by the value of this probability density function at its given coordinate. Let " & # * =! be the observed spatial coordinates in the data of the cells that exhibit the target phenotype. Then, ( ) is estimated via, where

is a Gaussian kernel with bandwidth parameter /, We then transform the kernel density estimates to fall between zero and one. Specifically, for cell , its color intensity would be computed via 1 # $%" exp 2 intensityH &I 6= 3H &I 2 min 3( $)

$*[)] max 3 ( $) 2 min 3 ( $) $*[)] $*[)] where

is the total number of cells in the dataset and [ ] are the set of integers from 1 to . To speed up computation, Monkeybread provides the option to run grid-based approximation 40. Specifically, a grid is overlayed onto the spatial data and each cell is assigned to the grid space that it falls within. Let ", & , + be the partitioning of cells into grid spaces and let , !, & , , " * =! be the coordinates of their centroids. Because cells within the same grid share the same <collapsed= coordinates, the kernel density estimate for a given location requires iterating over all grid spaces rather than all cells: 1 + $%" 3 ( ) 6= 7 | $| H 2 , #, /I In the limit, as the grid becomes fine-grained enough such that each grid space contains a maximum of one cell, this approximation will converge on the true kernel density estimate. A further speedup is made by coloring each cell according to the kernel density estimate at its grid space9s centroid rather than its own coordinate 3 htat is, the density estimate is calculated at each grid9s centroid and mapped back to each cell withinthat grid space. Lastly, to derive the color intensity used in the visualization, we transform these density values so that they fall between zero and one. Specifically, let ( ) be the function that maps each cell to the grid space that it falls within. Then, for cell , its color intensity is computed as intensityH &I 6= 3H -(&)I 2 min 3H -($)I

$*[)] max 3H -($)I 2 min 3H -($)I $*[)] $*[)] where

is the total number of cells in the dataset.

Monkeybread implements the statistical test for spatial co-expression described by He et al. (2021) 23. To review, this test considers the following test statistic: for a given ligand, , and receptor, , and two cell types,

and , under test, for each cell of type , we first identify the neighboring cells of type . We then compute the following ligand-receptor co-expression score: ( , ) 6=

V $,1 &,2 where of gene distribution.

is the set of all pairs of neighboring cells between types and , $,1 is the expression in cell (of cell type ), and &,2 is the expression of gene in cell (of cell type ).

Statistical significance of this score is assessed via a permutation test in which the pairings of neighboring cells are repeatedly permuted, and the test statistic is recalculated to form the null

The melanoma tumor dataset was downloaded from Vizgen9s website athttps://vizgen.com/datarelease-program and preprocessed with scanpy. Specifically, cells with less than 10 transcript counts were filtered out of the dataset resulting in a total of 189,071 cells. Transcript counts were normalized via log ( / + 1) where is the transcript count for a given gene in a given cell and are the total transcript counts in that cell. Before clustering, PCA was performed with 50 components and the k-nearest-neighbors (KNN) graph was computed on the principal components using k=15 neighbors. Cell types were annotated by first clustering the full set of cells using the Leiden clustering algorithm 22 on the KNN graph and then labeled according to their expression of known marker genes within each cluster (Figure S4). Clusters that expressed multiple marker genes (indicating potential heterogeneity with respect to their constituent cell types) were further sub-clustered and annotated (Figures S5-7).

Healthy Donors PBMNCs were acquired and CD16 negative monocytes (M0) were isolated with the Miltenyi biotec kit, (#130-117-337) following the manufacturer9s instructions. To develop polarized macrophages, monocytes were seeded on 100mm culture dishes with X-Vivo15 media supplied with 5% human serum. To generate M1 macrophages, cells were cultured for 5 days under GM-CSF stimulation and then for 2 days with IFNg and given a final 18h touch of LPS treatment. M2-like macrophages were developed through 5 days of culture with M-CSF and further 48h treatment of IL-4 stimulation. TAMdevelopment was accomplished via a continuous 7-day of cell culture under the influence of M-CSF+IL4 and IL10. The medium was enriched with 20% tissue conditioned medium (TCM) medium as the required supplement. Through the last 48h of culture, TGFb stimulation was added to the culture medium as well. Medium change for the culture systems was administrated once during the 7 days of culture. Tissue conditioned medium (TCM) was generated from MDA MB-231 breast cancer cell line following the protocol described by Benner et al. (2019) 32.

Cells derived from two patient donors were cultured in preparation for bulk RNA sequencing as described above, with four in vitro conditions per donor (monocyte control, M1 macrophages, M2 macropages, and TAM). After harvesting and counting, 1 million cells were aliquoted per culture condition and centrifuged at 400g for 5 minutes. Pelleted cells were resuspended in 350 uL Buffer RLT (Qiagen) + 1% BME, and vortexed for 30 seconds to 1 minute to disrupt the cells. Lysates were then frozen at -80# and stored until ready to proceed with RNA isolation. Total RNA was isolated from cell lysates using the RNeasy Mini Kit (Qiagen) per the <Purification of Total RNA from Animal Cells using Spin Technology= section of the kit manual. An optional DNase I treatment was included and performed using the RNase-Free DNase Set (Qiagen). Total RNA concentration was evaluated by RNA HS Qubit (ThermoFisher), followed by RNA ScreenTape (Agilent) to assess RNA quality. All RNA samples were of sufficient concentration and RNA Integrity (RIN) value to proceed with library construction.

Bulk libraries were prepared using the NEBNext Ultra ii Directional RNA Library Prep for Illumina Kit (New England BioLabs) per Section 1 of the manufacturer9s protocol. 500ng of total RNA was used as input for each sample, and an rRNA depletion was performed using the NEBNext Poly(A) mRNA Magnetic Isolation Module (New England BioLabs). SPRIselect beads (Beckman Coulter) were used in place of NEB Sample Purification Beads for all purification steps. For the final enrichment PCR, bulk libraries were barcoded using the NEBNext Multiplex Oligos for Illumina 3 96 Index Primers Kit (New England BioLabs) and amplified for 10 PCR cycles.

Bulk library quality was assessed by DNA HS Qubit assay (ThermoFisher) to measure library concentration, and by D5000 ScreenTape assay (Agilent) to evaluate library quality and size. All libraries were of sufficient concentration and quality to then prepare for sequencing. Library molar concentrations were more precisely determined by quantitative PCR (qPCR) using the KAPA Library Quantification Complete Kit 3 Universal (Roche). Based on calculated molar concentrations, bulk libraries were individually diluted to 20 nM and pooled equally by volume for sequencing.

Bulk libraries were submitted to Broad Institute Clinical Research Sequencing Platform for pairedend sequencing on an Illumina NextSeq500 system. The sequencing parameters consisted of 25 basepairs for Read 1, 50 basepairs for Read 2, and 8 basepairs for Index 1, with a target sequencing depth per sample of at least 20 million reads. Following sequencing, all raw data files were transferred to Immunitas servers for data analysis.

The preprocessing steps were carried out using a Snakemake pipeline 41,42. Briefly, the raw FASTQ files were aligned using STAR (v2.7.10b) 43 to genome GRCh38.primary_assembly. genome.fa using annotation file gencode.v40.annotaiton.gtf and the count matrix was generated by FeatureCounts (v2.0.1) 44. The raw counts were imported into DESeq2 45 (v1.36.0) in R (v4.2.1) and differential expression analysis was carried out following the standard procedure. Two differential expression analyses were performed to compute genes that were differentially expressed between M1-derived and M2-derived macrophages as well as between M1-derived and TAM-derived macrophages.

To generate cellular niches in the human melanoma tumor, we considered all immune cells and computed their normalized cell type neighborhood profiles using Monkeybread9s <neighborhoods_profile= function. All non-immune cells were labeled as <malignant/other=. We then computed the KNN graph using 100 nearest neighbors followed by clustering, via Leiden, with a resolution of 0.25. All clusters containing less than 300 cells were labeled as <malignant/other= whereas all remaining clusters were used as cellular niches for downstream analysis. This full pipeline is implemented in Monkeybread9s <cellular_niches= function.