Alexis J Combes
email: alexis.combes@ucsf.edu.

Bushra Samad

Jessica Tsui

Nayvin W Chew

Peter Yan

Gabriella C Reeder

Divyashree Kushnoor

Alan Shen

Brittany Davidson

Andrea J Barczak

Michael Adkisson

Austin Edwards

Mohammad Naser

Kevin C Barry

Tristan Courau

Taymour Hammoudi

Rafael J Argüello

Arjun Arkal Rao

Adam B Olshen

The Immunoprofiler consortium Cathy Cai

Jenny Zhan

Katelyn C Davis

Robin K Kelley

Jocelyn S Chapman

Chloe E Atreya

Amar Patel

Adil I Daud

Patrick Ha

Aaron A Diaz

Johannes R Kratz

Eric A Collisson

Gabriela K Fragiadakis

David J Erle

Alexandre Boissonnas

Saurabh Asthana

Vincent Chan

Matthew F Krummel
email: matthew.krummel@ucsf.edu

Discovering dominant tumor immune archetypes in a pan-cancer census

Keywords:

immunotherapy treatment strategies [START_REF] Binnewies | Understanding the tumor immune microenvironment (TIME) for effective therapy[END_REF][START_REF] Gajewski | Cancer immunotherapy strategies based on overcoming barriers within the tumor microenvironment[END_REF][START_REF] Gotwals | Prospects for combining targeted and conventional cancer therapy with immunotherapy[END_REF][START_REF] Hegde | The Where, the When, and the How of Immune Monitoring for Cancer Immunotherapies in the Era of Checkpoint Inhibition[END_REF].

To date, the combination of The Cancer Genome Atlas (TCGA) (Cancer Genome Atlas Network, 2015) and deconvolution technique base on gene expression from bulk tissue such as CIBERSORT [START_REF] Aran | xCell: digitally portraying the tissue cellular heterogeneity landscape[END_REF][START_REF] Chen | Profiling Tumor Infiltrating Immune Cells with CIBERSORT[END_REF][START_REF] Newman | Robust enumeration of cell subsets from tissue expression profiles[END_REF] have established a foundational but low resolution landscape of the TME across human tumors [START_REF] Bindea | Spatiotemporal dynamics of intratumoral immune cells reveal the immune landscape in human cancer[END_REF][START_REF] Gentles | The prognostic landscape of genes and infiltrating immune cells across human cancers[END_REF][START_REF] Mlecnik | Integrative Analyses of Colorectal Cancer Show Immunoscore Is a Stronger Predictor of Patient Survival Than Microsatellite Instability[END_REF][START_REF] Rooney | Molecular and genetic properties of tumors associated with local immune cytolytic activity[END_REF][START_REF] Thorsson | The Immune Landscape of Cancer[END_REF]. More recently, single-cell RNA sequencing (scRNA-seq) technologies have been increasingly applied to define the diversity of TME in a single cancer type [START_REF] Azizi | Single-Cell Map of Diverse Immune Phenotypes in the Breast Tumor Microenvironment[END_REF][START_REF] Goswami | Immune profiling of human tumors identifies CD73 as a combinatorial target in glioblastoma[END_REF][START_REF] Lavin | Innate Immune Landscape in Early Lung Adenocarcinoma by Paired Single-Cell Analyses[END_REF][START_REF] Zhang | Landscape and Dynamics of Single Immune Cells in Hepatocellular Carcinoma[END_REF] or focus on a single immune cell compartment [START_REF] Cheng | A pan-cancer single-cell transcriptional atlas of tumor infiltrating myeloid cells[END_REF][START_REF] Gueguen | Contribution of resident and circulating precursors to tumor-infiltrating CD8+ T cell populations in lung cancer[END_REF][START_REF] Oh | Intratumoral CD4+ T Cells Mediate Anti-tumor Cytotoxicity in Human Bladder Cancer[END_REF][START_REF] Zhang | Single-Cell Analyses Inform Mechanisms of Myeloid-Targeted Therapies in Colon Cancer[END_REF][START_REF] Zilionis | Single-Cell Transcriptomics of Human and Mouse Lung Cancers Reveals Conserved Myeloid Populations across Individuals and Species[END_REF]. The identification of recurrent motifs in the immune system at a higher resolution and spanning multiple cancer types is still lacking.

It is now well understood that complex coordination of immune cell states is required to achieve important tissue functions such as wound healing, tissue homeostasis or viral clearance [START_REF] Mujal | Immunity as a continuum of archetypes[END_REF]. Rather than focusing on the state of a single cell type, a given immune response can be conceived as a collection of cell subsets and specific immune cell pairings linked with function allowing us to define "immune archetypes". For response to therapy, a strong archetypal relationship between immune cell subsets such as cDC1, CD8 T cells and NK cells [START_REF] Barry | A natural killer-dendritic cell axis defines checkpoint therapyresponsive tumor microenvironments[END_REF][START_REF] Böttcher | NK Cells Stimulate Recruitment of cDC1 into the Tumor Microenvironment Promoting Cancer Immune Control[END_REF] or cDC2 and CD4 conventional and regulatory T cells [START_REF] Binnewies | Unleashing Type-2 Dendritic Cells to Drive Protective Antitumor CD4+ T Cell Immunity[END_REF][START_REF] Bosteels | Inflammatory Type 2 cDCs Acquire Features of cDC1s and Macrophages to Orchestrate Immunity to Respiratory Virus Infection[END_REF] have been previously identified in specific tumor types and viral infection. Tumors have also been described as 'Hot' (immune infiltrated) or 'Cold' (sparsely infiltrated) with variations in stromal content [START_REF] Bagaev | Conserved pan-cancer microenvironment subtypes predict response to immunotherapy[END_REF]. However, these are likely incomplete delineations and whether they are associated with tumor biology and granular immune cell composition is still unclear [START_REF] Chen | Elements of cancer immunity and the cancer-immune set point[END_REF].

In this study, we leveraged a unique dataset composed of both cell type compositional and transcriptomic data from 364 fresh surgical specimens across 12 tumor types to identify conserved tumor immune archetypes. We used an unsupervised clustering approach based on tumor specific immune gene signatures, benchmarked against 'ground truth' cell type composition data from flow cytometry, to identify and validate 12 unique tumor immune archetypes that were also identified in the TCGA dataset. These archetypes, discovered with only ten features, differentially aggregated other cell types that were not part of the discovery features, and delineate immune and tumor transcriptomic programs across different tissue type. They and these data thus provide an unprecedented resource to study cancer immunity and cancer targets to improve response to immunotherapy.

RESULTS

A Pan-Cancer High-Dimensional Study of Dominant Immune Composition.

The UCSF Immunoprofiler Initiative (IPI) collected fresh surgical specimens from 12 distinct tissues of origin using an unbiased approach, i.e., agnostic of tumor type, stage, and grade (Figure 1A, Supplementary Table 1). We performed standardized processing of 364 individual tumor specimens including rapid digestion into single cell suspension for immune phenotyping using multi-parametric flow cytometry (flow panels, see STAR methods). To identify patterns of gene expression within six broadly defined cell populations, we also performed cell sorting for bulk RNA-sequencing including compartments denoted: 1. "Live": All viable cells at the time of sorting, 2. "T conv ": sorted conventional CD4 + and CD8 + T cells, 3. "T reg ": CD25 + CD4 + (enriched for regulatory) T cells, 4. "Myeloid": Lymphocyte-negative HLA-DR + (enriched for myeloid) cells, 5. "Stromal": CD45 -CD44 + Thy1 + cells, and 6. "Tumor": all other CD45 -cells (Figure 1B and S1A left, and see STAR methods for full sort descriptions and potential caveats).This approach was chosen over single-cell RNA sequencing because greater read depths capture weaker transcripts and because this technology was an immature technology at the start of the program.

Because they are the basis for many previous single-cancer type studies, we initially focused upon the abundance of three major cell types in the TME by flow cytometry (hereafter "3-feature" or "3-f"), namely total α/β T cells (Tcell feature), myeloid cells (Myeloid feature) and non-immune CD45 -, CD44 + , Thy1 + stromal cells (CD44 + CD90 + Stroma feature) (Figure S1A left and supplementary table 2). Consistent with previous descriptions [START_REF] Bagaev | Conserved pan-cancer microenvironment subtypes predict response to immunotherapy[END_REF][START_REF] Mandal | The head and neck cancer immune landscape and its immunotherapeutic implications[END_REF][START_REF] Thorsson | The Immune Landscape of Cancer[END_REF][START_REF] Varn | Systematic Pan-Cancer Analysis Reveals Immune Cell Interactions in the Tumor Microenvironment[END_REF] the cell abundances for these three cell types, plotted as a score (Figure 1B, and extended method) for each tumor type show tremendous heterogeneity across and within tumors from different tissues, suggesting the need for a tumor classification that goes beyond tissue of origin (Figure 1C).

Taking advantage of a distinguishing feature of our cohort which have linked compositional flow cytometry and RNAseq data, we discovered and then validated a set of gene signatures that would allow us to infer and compare compositional data from external RNAseq-only datasets such as TCGA. We used differential gene expression (DGE) analysis on high quality RNA seq data from sorted tumor associated T cell, myeloid and stromal cells characterized by well separated clusters of transcriptomic profiles revealed by K-means clustering . We defined unique gene signatures for the Tcell, Myeloid and CD44+ CD90+ Stroma features, comprised of 25, 29 and 21 genes respectively, and used this to generate a compositional score where a high value corresponded with high cell abundancy in the tumor specimen analyzed (Figure 1D and Figure S1B ). Use of this score, when applied to RNAseq data from the live RNAseq compartment in a validation cohort showed high concordance with the cell type frequencies obtained via flow cytometry, independent of tissue of origin (Spearman correlation 0.91,0.90 and 0.94 of T cell, Myeloid, and Stromal, respectively) (Figure 1D). Moreover, rank ordering of tumor from different tissues by these gene scores showed similar trends, and strong correlations across individuals to those obtained by flow cytometric measurements (Figure 1C/S1D and 1E).

Our signatures for intra-tumoral composition significantly outperformed signatures derived from cell lines ( e.g. CIBERSORT Figure S1C) when compared to the 'ground-truth' of flow cytometric composition data. We presume this improvement results from our data being taken from cells directly collected from fresh tumor tissue as opposed to cell lines or isolated samples [START_REF] Aran | xCell: digitally portraying the tissue cellular heterogeneity landscape[END_REF][START_REF] Newman | Robust enumeration of cell subsets from tissue expression profiles[END_REF]. When we applied our gene signatures to RNA-seq data from 4341 TCGA tumor specimens RNA-seq data, we found that the median score for each tumor type between datasets for the majority of the tissues surveyed was similar, suggesting that the abundancies we describe by tumor type ( e.g. Figure 1C) extend beyond our sample processing protocol (Figure S1E). The exceptions to this were for lung (LUNG), liver (HEP), Glioblastoma (GBM) and pancreatic (PDAC) tumors (Figure S1E), this may be due to incomplete or different methods of sampling or patient-selection variation between the cohorts; thus, in the remainder of this analysis, we did not consider these cancer type when making additional comparisons.

Unsupervised Clustering of 3-Features, Independent of Tissue of Origin

We performed unsupervised clustering, using the Louvain community detection algorithm on a K nearest-neighbor (KNN) weighted graph, on 3 features, Tcell, Myeloid and CD90+ CD44+ Stromal scores for all 260 samples with a "Live" bulk RNAseq compartment (Supplementary Table S2). The clustering was visualized using UMAP, both on the IPI and TCGA cohorts (Figure 2A/E and STAR methods). The optimal clustering parameters were evaluated by minimizing the Davies Bouldin Index (DBI) (Davies and Bouldin, 1979), a metric that assesses the ratio of the intra-cluster distance to inter-cluster distance (Figure S2A-B and STAR methods) and in both cohorts this resulted in six clusters. In the IPI cohort, this separated two immune rich clusters (termed immune rich (IR) and immune stromal rich (IS)) defined by high expression of Tcell and Myeloid features and differentiated by Stromal enrichment in IS. This similarly differentiated two Immune desert (ID) clusters (Immune Desert and Immune Stromal Desert) defined by low expression of Tcell and Myeloid features, and again differentiated by enrichment of CD44 + CD90 + stroma within the Immune Desert cluster. Finally, we also identified clusters with enrichment in only one immune feature, namely the Tcell Centric (TC) and Myeloid Centric (MC) features . Application of our gene scores to the TCGA cohort also produced six clusters when minimizing the DBI and these cluster corresponded to the six clusters found in the IPI cohort (Figure 2E-F and Figure S2C-D).

To provide validation of the flow cytometry data, we ran immunofluorescence assays on tissues collected from a subset of the IPI tumor surgical specimens. Consistent with expectation, we observed that samples taken from immune rich clusters (1-3) had the highest CD45 + infiltration, while immune desert samples had the lowest (Figure 2G-H and Figure S2G-L). Myeloid Centric (MC) clusters were similar to the immune desert clusters for overall CD45+ infiltration as assessed by immunofluorescence (Figure 2G-H and Figure S2G-L). However, MC clusters presented an intermediate CD45+ cells frequency when measured by flow cytometry. Overall, we observed significant correlation between the frequency of immune cells among all cells versus among viable cells after tissue dissociation (Figure S2F). The discrepancy between flow cytometry and immunofluorescence observed in MC archetype may represent a modest preferential recovery of immune cells during tissue dissociation for MC tumors. Future analyses of spatial dimensions may reveal paired infiltration patterns i.e. T cell plus or minus myeloid or 'excluded' versus 'infiltrated' [START_REF] Galon | Type, density, and location of immune cells within human colorectal tumors predict clinical outcome[END_REF][START_REF] Mlecnik | Integrative Analyses of Colorectal Cancer Show Immunoscore Is a Stronger Predictor of Patient Survival Than Microsatellite Instability[END_REF]. This will involve considerable orthogonal analyses due to the intrinsic complexity of diverse tissue morphologies e.g. lung versus colon.

Given the importance of chemokines in recruitment of immune cells, we sought to determine whether these archetypes had unique chemokine gene-expression signatures that might corroborate their classification. We first derived a single gene signature score applied to the live compartment, based on 39 chemokines [START_REF] Nagarsheth | Chemokines in the cancer microenvironment and their relevance in cancer immunotherapy[END_REF], and found a high score in immune rich clusters and a low score in immune desert cluster in both the IPI and TCGA cohorts (Figure 2K-L and Figure S2E). Two exceptions stood out using the pan-chemokine measure: Tcell and Myeloid Centric clusters showed decreased and increased pan-chemokine scores, respectively (Figure 2K-L). We thus examined chemokines individually, using hierarchical clustering of the median chemokine expression amongst the tumors of each archetype, and identified sets of chemokines associated with each cluster (Figure 2M-N). For instance, the three T cell enriched clusters are enriched in either chemokines expressed by T cell (XCL1, XCL2, CCL4) or T cell attracting chemokines (CXCL13, CCL5). Other archetypes had their own pattern of differentially expressed chemokines and these were broadly consistent with their composition. Immune desert clusters typically had a few unique and specific chemokines expressed. The overall patterns observed on a chemokine-by-chemokine basis for the six archetypes was similar between IPI and TCGA (Figure 2N), with a few exceptions such as CCL1, CCL11 and CCL27. These results demonstrate that a 3-feature scoring of tumor tissue identifies six unique clusters defined by a different degree of immune infiltration and expression of distinct sets of chemokines.

Distribution of 3-feature Archetypes by Cancer Type and outcome

3-feature based clustering demonstrated a heterogenous distribution of tumor types among the different clusters, although some cancer type, such as Kidney and Melanoma had significant biases (Figure 3A and Figure S3A). Taking advantage of the large clinical dataset available in the TCGA cohort, we sought to broadly assess whether these simple archetypal classifications had a relationship to prognosis. Agnostic to the tissue of origin, the overall survival at 5 years analyzed by multivariate regression is significantly better for the immune rich archetype compared to all the others (p-value 3.9E -11 ) with a general trend in outcome that tracks with overall immune infiltration (Figure 3B). To extend this analysis and focus within tumor types, we analyzed distributions of archetypes in both IPI and TCGA and assessed the outcome in individual cancer (Figure 3C and Figure S3B-D). While the relative composition of archetypes by tissue of origin were broadly similar in both cohorts, there was some variation e.g. heavily 'Tcell Centric' bias in Melanoma in IPI cohort can appear as a combination of 'IR' or 'IS' in TCGA. Furthermore, archetypical classification in some tumor types such as melanoma, kidney, bladder, and colorectal appear to stratify survival outcome, trends were weaker within other individual cancer type. This variability prompted us to consider that these coarse-grained 3-feature archetypes may not capture the full range of immune cell heterogeneity in tumors. Measurement by flow cytometry of Treg density and CD4 + /CD8 + conventional T cell ratio confirmed our assumption as we observed variability within each archetype (Figure 3D-E).

Developing A 6-Feature Archetype Definition

Using flow cytometry data, we again grouped tumors by cancer type to assess variation in regulatory T cells (Treg) in live cells and CD4 + or CD8 + conventional T cell frequencies in T cells, both between and within tumor types (Figure 4A) (supplementary table 3). Repeating our previous methodology (Figure S1), we used DGE between the Treg RNAseq compartment and the other cell sorted RNAseq compartments to generate Treg gene signature composed of 9 genes, (Figure S4A) that was highly similar to other published signatures [START_REF] Vargas | Fc Effector Function Contributes to the Activity of Human Anti-CTLA-4 Antibodies[END_REF][START_REF] Plitas | Regulatory T Cells Exhibit Distinct Features in Human Breast Cancer[END_REF][START_REF] Zemmour | Single-cell gene expression reveals a landscape of regulatory T cell phenotypes shaped by the TCR[END_REF]. To isolate a signature for CD8 versus CD4 within our T conv RNAseq compartment, we used the flow data to identify samples rich in CD4 or CD8 conventional T cells and performed DGE between them (Figure S4B and STAR methods). Notably, most of the identified genes have been previously associated with CD4 + /CD8+ conventional T cell identity (e.g., CD8A, IL7R) or tissue residency (e.g., VCAM1, BACH2) [START_REF] Richer | T Cell Fates Zipped Up: How the Bach2 Basic Leucine Zipper Transcriptional Repressor Directs T Cell Differentiation and Function[END_REF][START_REF] Savas | Single-cell profiling of breast cancer T cells reveals a tissueresident memory subset associated with improved prognosis[END_REF] (Figure S4B). Assessment of a validation cohort (samples that were not used to generate the features scores) showed very high correlations (Spearman correlation Treg :0.86, CD4:0.97 and CD8:0.98, respectively), with population abundances measured by flow cytometry. Again, these correlations were significantly better than those obtained using CIBERSORT (Figure S4C). Unsupervised clustering of all samples using the CD4 and CD8 feature signature genes within our T conv RNAseq revealed the existence of at least 3 distinct groups of tumors across cancer type: CD8-biased tumors, CD4-biased tumors and a large population that was mixed for both these signatures (Figure 4B).

Using these scores, we next performed 6-feature clustering (Figure 4C and Figure S4D-F). DBI optimization yielded eight clusters and analysis by alluvial plot revealed that the two new clusters formed are mostly a subdivision of the previous coarse Immune Rich and Immune Stromal Rich archetypes, now delineated by CD4 to CD8 ratio (Figure 4D-E and Figure S4E). This increase in feature granularity also resulted in some specimens shifting within the classification. For example, some Tcell Centric samples now shifted to being considered "Immune-rich:CD4" because of their profound high CD4 score, Furthermore, this analysis revealed that the 3-feature cluster of 'Immune Stromal Desert' was relatively CD8 rich whereas the 'Immune Desert' favored CD4 cells. CD45 densities remain high for the immune rich tumors (Figure S4G). Assessing chemokine gene expression showed that this re-clustering dramatically segregated chemokine gene expression in the Immune Rich archetype and further refined the chemokine expression found amongst the other archetypes (Figure 4F and Figure S4H)

Mapping T Cell Exhaustion and Myeloid Subset Heterogeneity In 6-Feature Archetypes

Exhaustion in T cells (T ex ) represents a transcriptional state for T cells that arises in cancers and chronic viral infections and is characterized by progressive loss of effector functions, high and sustained inhibitory receptor expression, and acquisition of a distinct transcriptional program [START_REF] Blank | Defining "T cell exhaustion[END_REF]. We used the T conv RNA-seq compartment to identify genes that showed the highest correlation with CTLA4, PDCD1, HAVCR2, CD38 and LAG3, previously identified exhaustion markers. We identified 11 such genes, which included TOX, a transcription factor recently described as key driver of T cell exhaustion [START_REF] Beltra | Developmental Relationships of Four Exhausted CD8+ T Cell Subsets Reveals Underlying Transcriptional and Epigenetic Landscape Control Mechanisms[END_REF][START_REF] Khan | TOX transcriptionally and epigenetically programs CD8+ T cell exhaustion[END_REF]Scott et al., 2019) (Figure S4I and STAR methods). We then benchmarked this gene signature score against the abundance of CD4 + , CD8 + and the sum of both CD4+ and CD8+ T cells co expressing CTLA4, PD-1 and CD38 as markers of T ex (Figure S4J). This T ex gene score best mirrored the flow-cytometry-based abundance of CTLA-4+/PD-1+ cells within the combined CD8+ and CD4+ compartments (Figure S4K). Using this T ex gene signature score we observed enrichment of exhaustion in CD8 biased archetypes including the immune desert (ID) archetype (Figure S4L-M) and together with low MHC I expression by the tumor cells in these archetypes (Figure S4N).

In addition to T cell exhaustion, we also sought to assess myeloid heterogeneity in our tumor landscape, widely known to be variable in tissues. We thus probed the abundance of mononuclear phagocytes (MNP) subsets including monocytes (Mo), macrophages (Mp), classical dendritic cells (cDC2, cDC1) and plasmacytoid dendritic cells (pDCs) for each archetype using flow cytometry (Figure S4O). This revealed that, despite, slight enrichment between archetypes e.g., of cDC2 over cDC1 in IR CD4 Bias, the cDC1/cDC2 ratio, monocyte, and macrophage densities were highly variable within each of the 6-feature archetypes. (Figure 4G and Figure S4P).

10-Feature Archetype Definition

The complexity of the MNP identity and plasticity has complicated efforts to determine which populations are beneficial or subversive to the anti-tumor response [START_REF] Broz | Dissecting the tumor myeloid compartment reveals rare activating antigen-presenting cells critical for T cell immunity[END_REF][START_REF] Etzerodt | Tissue-resident macrophages in omentum promote metastatic spread of ovarian cancer[END_REF][START_REF] Gubin | High-Dimensional Analysis Delineates Myeloid and Lymphoid Compartment Remodeling during Successful Immune-Checkpoint Cancer Therapy[END_REF]. To assess myeloid diversity and density across all samples, we generated a scRNAseq sub-study of sorted tumor associated MNP to identify specific gene signatures for the principal MNP subsets (Figure S5A). Following removal of non-MNP cellular contaminants, we found 5 unique clusters from 3,880 input cells (Figure S5B). Using DGE between each cluster, we identified unique 5 gene signatures for each subset (Figure S5C). Most of the genes identified have been previously associated with their respective MNP subsets such as VCAN, TREM2, CD1C, CLEC9A and LAMP5 in monocytes, macrophages, cDC2, cDC1 and pDCs, respectively [START_REF] Binnewies | Understanding the tumor immune microenvironment (TIME) for effective therapy[END_REF][START_REF] Cheng | A pan-cancer single-cell transcriptional atlas of tumor infiltrating myeloid cells[END_REF][START_REF] Combes | BAD-LAMP controls TLR9 trafficking and signalling in human plasmacytoid dendritic cells[END_REF][START_REF] Molgora | TREM2 Modulation Remodels the Tumor Myeloid Landscape Enhancing Anti-PD-1 Immunotherapy[END_REF][START_REF] Sancho | Identification of a dendritic cell receptor that couples sensing of necrosis to immunity[END_REF]. To validate these signatures across individuals and tissue types we performed a similar analysis on 11 fresh resected tumors across 3 tissue types, Melanoma, Head and Neck and Kidney (Figure S5D). We confirmed that the same 5 gene signatures, that we observed in our initial melanoma sample, uniquely identified the MNP subsets independent of the tumor type (Figure S5D-E). A sixth cluster characterized by high expression of FSCN1, CCR7 and LAMP3 was also detected in this dataset (Figure S5E). This specific gene expression signature has been recently associated with a distinct cDC transcriptomic/activation state (provisionally named 'cDC3') characterized by high expression of regulatory molecules such as PD-L1 and highly correlated with Treg abundance in the tumor [START_REF] Cheng | A pan-cancer single-cell transcriptional atlas of tumor infiltrating myeloid cells[END_REF][START_REF] Maier | A conserved dendritic-cell regulatory program limits antitumour immunity[END_REF][START_REF] Mulder | Cross-tissue single-cell landscape of human monocytes and macrophages in health and disease[END_REF][START_REF] Zhang | Landscape and Dynamics of Single Immune Cells in Hepatocellular Carcinoma[END_REF]. However, due to the absence of specific markers in our antibody panel to validate the presence of this subset by flow cytometry we could not include this population subset in our MNP features scores.

Gene signature scores for these MNP subsets were generated in bulk RNAseq for tumor associated myeloid cells and validated against cellular abundance by flow cytometry. The correlation of each MNP subset feature score across all cancer type were highly correlated (Spearman correlation of 0.86, 0.91, 0.93, 0.94 and 0.94 for Monocytes, Macrophages, cDC2, cDC1 and pDC, respectively) (Figure S5F and STAR methods). Again, when we analyzed the distributions of these five MNP subset features scores by cancer type, a high heterogeneity was apparent (Figure 5A and Figure S5G). Thus, we repeated the unsupervised clustering after adding four MNP subset features including macrophages, monocytes and the two types of classical dendritic cells to the previous six features (Supplementary Table 3). The resulting clustering (hereafter "10feature" or "10-f") produced 12 unique tumor immune archetypes after DBI minimization, with no bias toward specific archetypes found between scores calculated from flow and scores calculated directly from RNAseq (Figure 5A and Figure S5H-J). As previously observed when adding T cell subsets, the addition of the four MNP measures subdivided preexisting immune archetypes but did not increase those by 16-fold as could occur if these were randomly assorted (Figure 5B-C). MNP inclusion also resulted in some specimens shifting between MC and ID archetypes, driven by strong monocytes or macrophages enrichment.

Analyzing the 12 archetypes produced in 10-f clustering, we found that generally both IR and ID sub-archetypes were distinguished by distinct pairings of T cell and myeloid subsets but that these pairings varied. Specifically, the monocyte/macrophage ratio demarcates two different IR CD8 bias archetypes where Treg abundance is generally higher in macrophageenriched tumor specimen (Figure 5D-F and 5O-Q). One archetype, IR CD4 biased, is differentiated from its IR CD8 counterparts by enrichment in cDC2 vs cDC1 (Figure 5D and 5R-T). However, such correlation between CD4/CD8 and cDC2/cDC1 ratios is only observed among IR archetypes and not in archetypes containing high stromal densities (Figure 5G-I and 5R-T).

The relationship between monocytes/macrophages and T cell subsets composition also differed significantly amongst the 12 archetypes (Figure 5D-N). In IS (Figure 5H), CD8 abundance is opposed to Macrophage abundance, but in ID, both CD4 and CD8 rich sub-archetypes are equivalently enriched in macrophages and monocytes (Figure 5N, 5P and 5T). TC (Tcell centric) archetypes are characterized by a positive correlation between Treg and Macrophage abundance . This discordance between the abundances of macrophages and other immune cells in different archetypes could suggest that a generalized Macrophage score does not capture critical heterogeneity in macrophage phenotypes or that macrophages and T cells may engage in additional as-yet-undiscovered relationships that regulate their numbers.

Defining Immune Gene Expression Pattern Across 10-Feature Archetypes

Focusing on the 12 tumor archetypes identified by 10-f clustering, we examined associations with cell abundances and gene sets that were not used for clustering. Significant archetypespecific enrichment was found for gene sets that define intra-tumor NK cells in IR based on published signatures [START_REF] Barry | A natural killer-dendritic cell axis defines checkpoint therapyresponsive tumor microenvironments[END_REF] and these were confirmed by flow cytometry measurement, while ISR were enriched in gene sets defining Mast cells [START_REF] Cheng | A pan-cancer single-cell transcriptional atlas of tumor infiltrating myeloid cells[END_REF] (Figure 6A and S6B-C). Plasma cells, defined by enrichments of IgG genes, were enriched in T cell Centric Macrophage bias tumors whereas B cell measured by proteomic and transcriptomic expression (MS4A1; CD20) were most prevalent in T cell Centric DC rich [START_REF] Chen | Single-cell transcriptome and antigen-immunoglobin analysis reveals the diversity of B cells in non-small cell lung cancer[END_REF]van Galen et al., 2019) (Figure 6A and S6D). These may correspond to a prevalence for tertiary lymph nodes in those archetypes although we did not observe higher DC frequencies in those tumor specimens (Figure S5J). An archetype-specific immune composition pattern again correlated with specific enrichment of groups of chemokine transcripts (Figure 6B). Notably, chemokines with specificity for families of chemokine receptors now frequently co-clustered within the 10-f archetypes.

We next used the myeloid RNAseq compartment to explore the level of expression of genes previously associated to tumor-infiltrating MNP, namely 'M1', 'M2', 'DC' and 'Costimulatory molecules' [START_REF] Biswas | Tumor-associated macrophages: functional diversity, clinical significance, and open questions[END_REF][START_REF] Cassetta | Human Tumor-Associated Macrophage and Monocyte Transcriptional Landscapes Reveal Cancer-Specific Reprogramming, Biomarkers, and Therapeutic Targets[END_REF][START_REF] Cheng | A pan-cancer single-cell transcriptional atlas of tumor infiltrating myeloid cells[END_REF][START_REF] Maier | A conserved dendritic-cell regulatory program limits antitumour immunity[END_REF][START_REF] Roberts | Critical Role for CD103(+)/CD141(+) Dendritic Cells Bearing CCR7 for Tumor Antigen Trafficking and Priming of T Cell Immunity in Melanoma[END_REF]. Despite the heterogeneity in individual gene expression associated to 'M1' and 'M2' macrophage phenotype among archetypes we observed a general enrichment in 'M2' genes in IS archetypes (Figure 6C). On the other hand, 'DC' and 'Co-Stim' associated genes were highly expressed in the Tcell Centric DC rich archetype but LAMP3, FCSN1 and CCR7 expression (typically indicative of the "cDC3" state) did not appear coordinated across archetypes (Figure 6C). This may indicate a heterogeneity of dendritic cell transcriptomic states between archetypes.

To further elucidate the heterogeneity found in myeloid function we combined the myeloid 'signature' gene sets with a selection of 12 other gene sets highly linked to subsets of Treg, T conv , B cell, plasma cells, and NK cells and performed hierarchical clustering of these genes sets based on their median expression in the archetypes. The gene sets were evaluated in different RNAseq compartments corresponding to their appropriate cell type, namely 'Th1' genes in the T conv RNAseq compartment, 'M1' genes in the myeloid RNAseq compartment and 'Treg homing' in the Treg RNAseq compartment (Figure 6D and Supplementary table 3). Combining the gene sets in this way revealed distinct immune signatures for each archetype (Figure 6D). For instance, while the IR CD8 macrophage bias archetype was enriched in genes associated with the type 1 response (IFNG, TNF, IL1B) both in T conv and MNP compartments, both IS and ID CD8 archetypes were characterized by their unique combination of upregulated gene expression associated to T cell exhaustion (PDCD1, CTLA4, ENTPD1) together with 'M2' like macrophages (PDCD1L2, CD163, CD274, and MRC1). This analysis also identified putative gene expression interactions present in the same archetype across cell types, such as high CCR2 expression by Tregs in an archetype rich in monocytes, which are known producers of the CCR2 ligand, CCL2 [START_REF] Loyher | Macrophages of distinct origins contribute to tumor development in the lung[END_REF][START_REF] Mondini | CCR2-Dependent Recruitment of Tregs and Monocytes Following Radiotherapy Is Associated with TNFα-Mediated Resistance[END_REF] (Figure 6D). Taken together, this demonstrated that the 12 tumor immune archetypes found in 10-f analysis also identifies a unique combination of cell composition and transcriptomic profiles.

Notably, certain immune populations remained more variable across the archetypes. For instance, neutrophils were mostly uncorrelated although there was a slight and statistically insignificant rise in neutrophil infiltration in archetypes with low immune infiltration (Figure S6A and S6E). A gene signature previously associated to stimulatory dendritic cells abundance and better outcome in Melanoma patients is slightly enriched in archetypes with high total immune infiltration [START_REF] Barry | A natural killer-dendritic cell axis defines checkpoint therapyresponsive tumor microenvironments[END_REF]. Exhaustion signature is enriched in both immune rich and immune desert archetypes (Figure S6A and S6F-G). Thus, while the 12 immune archetypes from 10-f analysis identify dominant combinations of cell types and transcriptomic profiles present in tumors, other cell types may also be independently layered within and across this apparently conserved biology.

10-Feature Tumor Archetypes Tie Closely to Tumor Biology and Disease Outcome

Finally, we sought to examine the relationship of 10-f archetypes to the phenotype of the tumor cells themselves. As previously shown for the coarse archetypes, 10-f archetypes are diverse in their tissue of origin; however, some tumor types remain highly represented within an archetype (e.g., Kidney: Immune Rich CD8, Melanoma: T cell centric archetypes) (Figure 7A and Figure S7). Regardless of the frequent intra-cancer heterogeneity, we observed a strong relationship between tumor proliferative capacity, measured by the fraction of tumor cells which were Ki67 + as measured by flow, and immune desert and myeloid centric archetypes (Figure 7B).

Based on this we assembled gene signatures for various aspects of tumor biology to determine whether immune features mapped to tumor biology. Expression of cell cycle associated genes indicated enhanced tumor cell proliferative capacity in ID archetypes, consistent with the Ki67 data (Figure 7C). Conversely, IR and ISR tumors were enriched in other transcriptomic programs such as Interferon-Stimulated Genes ISG, Senescence or Epithelial-mesenchymal Transition (EMT) related genes [START_REF] Kinker | Pan-cancer single-cell RNA-seq identifies recurring programs of cellular heterogeneity[END_REF][START_REF] Muñoz | Targetable mechanisms driving immunoevasion of persistent senescent cells link chemotherapy-resistant cancer to aging[END_REF][START_REF] Wiley | Analysis of individual cells identifies cell-to-cell variability following induction of cellular senescence[END_REF] and only some immune deserts were highly enriched in fibrosis associated genes (Figure 7C). Without pre-selection of genes from previous literature, we performed DGE and then hierarchical clustering using the most differentially expressed genes, between archetypes on the sorted tumor RNAseq compartment (Figure 7D). Within the 12 sets of differential genes, we noted some discernable patterns, for example increased IFNG expression in tumor cells from the IR CD8 Macrophage archetype, concordant with increased type 1 response in immune cells from the same archetype.

To define gene signatures that we could assess in TCGA, we queried our sorted live compartment to find genes that correlated with each archetype (Figure 7E). This identified unique gene signatures for each of the 12 archetypes, within which we observed both immune (e.g. LAYN, CTLA4, CSF1, CCL19) and non-immune (CRHBP, PTGFR, HSPA2, MTCL1) related genes. Applying these signatures to the TCGA dataset, we found that the tumor archetypes identified in our IPI dataset were retrieved with a similar relative composition of archetypes by cancer type in KIRC, CRC, LUAD, HNSC, BLCA or with a slight shift between archetypes from the same coarse archetype (i.e. : ID CD4 Bias to IR CD8 bias) in GYN and SKCM (Figure S7A-K). Additionally, we observed a similar distribution of archetypes in both cohorts when looking at clinical correlates such as tumor stage, grade, and metastatic status (Figure S7L-N). This suggests that our archetypes can be identified in independent datasets and the predictive potential of these signatures for rapid classification of patients in both primary and metastatic tumors independent of tumor grade and stage.

Finally, survival analysis of the different tumor archetypes identified in the TCGA dataset for each cancer type separately showed that better outcomes may be cancer type specific, and each cancer type may have a different type of TME promoting the best immune response (Figure S7A-K). However, when survival analysis was performed across tumors using a multivariate survival regression, we detected significant outcome differences between archetypes that have similar T cell subset enrichment regardless of the tissue of origin (Figure 7F). For instance, in IR CD8 archetypes the apparent enrichment in monocytes over macrophages (Pink archetype versus Red) is associated to better survival (Median survival:3354 days) and more so when compared to ID and IS CD8 biased archetypes (p-value 1.55E-12) (Figure 7F top). Conversely, in archetypes with no significant bias between CD4 and CD8 T cell, TC archetypes characterized by an enrichment of Macrophages over monocytes displayed the higher median survival (2456 days, p-value 1.17E-6) while in archetypes with CD4 T cell bias no significant outcome difference between archetypes was detected regardless of myeloid biases (p-value 3.01E-1) (Figure 7F middle and bottom). This indicates that tumor archetypes provide a template to study different subtypes of anti-tumor immune response in a variety of primary human cancers. Finally, our study represents an initial framework-one certain to be refined-to understand how tumor archetypes relate to other biological contexts (Supplementary table 6) and how to best modulate these archetypes depending on their specific immune context.

DISCUSSION

In this study, we present a holistic survey of dominant immune archetypes across 12 cancers from different tissues using fresh surgical tumor specimen from 364 patients. Empowered by complementary profiling assays we were able to discover 12 distinct immune archetypes that span cancer types. Each archetype is made up of a unique combination of cell composition -some used to derive the cluster, but many that are learned from it-and immune and tumoral transcriptomic phenotypes.

Starting with just 10 independent cell compositional features, unsupervised clustering revealed only 12 distinct clusters whereas 10 independently assorting binary variables might have produced 2 10 or 4096. While more work will need to be done to ensure that this is not due to a lack of sampling, this implies that some combinations of cell densities may not exist in the TME of solid tumors. This is partly in line with previous work which identified between four and six immune subtypes spanning multiple tumor types in TCGA using mainly immune cell fractions from deconvolution analysis of all tissue bulk RNA-seq data and knowledge-based immune gene expression signatures [START_REF] Bagaev | Conserved pan-cancer microenvironment subtypes predict response to immunotherapy[END_REF][START_REF] Thorsson | The Immune Landscape of Cancer[END_REF]. Notably, similar dominant immune pathways are also revealed in our analysis such as the IFN pathway or the imbalance between T cell and MNP subsets. However, these similarities go beyond immune related features. We also identified stromal frequency as a key distinctive feature in the tumor classification across tissues that subdivides both immune rich and immune desert archetypes. Additional features may partly explain our higher number of archetypes, but it may also lie in the power of a dataset that has large numbers of members and spans a diverse spectrum of cancers. Nevertheless, our progressive analysis from three coarse features through more granular definitions of T cells and MNP subsets suggest that it is the intersection of these 10 features that seem to arrive at a stable view of the major archetypes. However, further refinement and discovery is required especially in immune deserts where the paucity of overall immune cells likely decreases our ability to resolve differences using only these 10 measurements.

As previously shown by others, in a single cancer type or in tumor mouse models, the TME can be coarsely categorized into immune rich and immune desert areas [START_REF] Bagaev | Conserved pan-cancer microenvironment subtypes predict response to immunotherapy[END_REF][START_REF] Duan | Turning Cold into Hot: Firing up the Tumor Microenvironment[END_REF][START_REF] Galon | Approaches to treat immune hot, altered and cold tumours with combination immunotherapies[END_REF][START_REF] Mariathasan | TGFβ attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells[END_REF]. Our analysis revealed that the tumor immune archetypes subdivide this broad classification and identify distinct immune networks for each archetype with unique relationships amongst cell densities and chemokine networks. While IR archetypes are characterized by patterns previously described in Melanoma, Head and Neck or Breast malignancies-such as a correlation between CD8+ T cells, cDC1 ratio and NK cell abundance [START_REF] Barry | A natural killer-dendritic cell axis defines checkpoint therapyresponsive tumor microenvironments[END_REF][START_REF] Salmon | Expansion and Activation of CD103(+) Dendritic Cell Progenitors at the Tumor Site Enhances Tumor Responses to Therapeutic PD-L1 and BRAF Inhibition[END_REF] or between CD4+ T cell and cDC2 [START_REF] Binnewies | Unleashing Type-2 Dendritic Cells to Drive Protective Antitumor CD4+ T Cell Immunity[END_REF][START_REF] Michea | Adjustment of dendritic cells to the breast-cancer microenvironment is subset specific[END_REF]-other archetypes implicate different immune cell networks such as the apparent enrichment of plasma cells and Tregs in TC archetype. That particular co-association was previously shown to support the plasma cell residency in bone marrow and the presence of both cell types together associated with better response to immune checkpoint blockade therapy in sarcoma (Glatman [START_REF] Zaretsky | T Regulatory Cells Support Plasma Cell Populations in the Bone Marrow[END_REF][START_REF] Petitprez | B cells are associated with survival and immunotherapy response in sarcoma[END_REF]. It is tempting to propose that tumors hijack specific immune archetypes that originate for use in a completely different setting (Supplementary table 6).

This idea is perhaps also suggested by the distinct patterns of chemokine expression identified for each tumor immune archetype. For example, while IR CD8 archetypes are defined by a CXCL9, CXCL10, CXCL11/CXCR3 axis which resembles the biology of an ongoing chronic viral infection [START_REF] Metzemaekers | Overview of the Mechanisms that May Contribute to the Non-Redundant Activities of Interferon-Inducible CXC Chemokine Receptor 3 Ligands[END_REF], ID CD8 Macrophage bias archetype is defined by a strong expression of chemokines binding CX3CR1, which is an axis essential in settings of immune surveillance and homeostasis [START_REF] Gerlach | The Chemokine Receptor CX3CR1 Defines Three Antigen-Experienced CD8 T Cell Subsets with Distinct Roles in Immune Surveillance and Homeostasis[END_REF]. Further work is certainly needed to profile the vast array of other immune niches in homeostasis and immune responses, to test how these linked cell states are conserved across tissue and disease; the increasing availability of multi-omics resources will be an important step forward [START_REF] Mulder | Cross-tissue single-cell landscape of human monocytes and macrophages in health and disease[END_REF][START_REF] Reynolds | Developmental cell programs are co-opted in inflammatory skin disease[END_REF].

Tumor-specific genetics and mutational burden have been proposed to be key for anti-tumor immunity in multiple cancers [START_REF] Brown | Neoantigens predicted by tumor genome meta-analysis correlate with increased patient survival[END_REF][START_REF] Ghosh | Mutant p53 suppresses innate immune signaling to promote tumorigenesis[END_REF][START_REF] Goodman | Tumor Mutational Burden as an Independent Predictor of Response to Immunotherapy in Diverse Cancers[END_REF]. However, further investigation on extrapolating immune gene signature across TCGA cohort showed no correlation between gene expression and mutational burden in any cancer type [START_REF] Bagaev | Conserved pan-cancer microenvironment subtypes predict response to immunotherapy[END_REF][START_REF] Spranger | Density of immunogenic antigens does not explain the presence or absence of the T-cell-inflamed tumor microenvironment in melanoma[END_REF]. Conversely, we demonstrated that immune tumor archetypes are associated with tumor proliferation, diverse tumor transcriptomic programs and overall survival. Studying the synergy between immune archetypes and tumor mutational profile, may further elucidate the importance of tumor mutational burden in anti-tumor immune response [START_REF] Berger | A Comprehensive Pan-Cancer Molecular Study of Gynecologic and Breast Cancers[END_REF][START_REF] Quigley | Genomic Hallmarks and Structural Variation in Metastatic Prostate Cancer[END_REF].

In summary, our comprehensive characterization of the TME across many human solid cancer types reveals the existence of common and reproducible immune archetypes defined by distinct cell networks. To facilitate usage of our data for the wide research community, we will provide access to transcriptomic and compositional data for each archetype (datalibrary.ucsf.edu/public-resources). It is still unclear whether improving tumor cure rates will be based on the enhancement of a specific archetype, and the expansion of this classification in metastatic tumors as well as in animal models may help to test this [START_REF] Maynard | Therapy-Induced Evolution of Human Lung Cancer Revealed by Single-Cell RNA Sequencing[END_REF]. Moreover, it will be important to expand this analysis to the peripheral immune system as well as biopsies before and after immunotherapies in order to define the relationship between 'dominant' and 'reactive' tumor immune archetype across tumor types [START_REF] Bi | Tumor and immune reprogramming during immunotherapy in advanced renal cell carcinoma[END_REF]Grünwald et al., 2021). Furthermore, there will undoubtedly be many additional ways to explore this dataset and discover patterns of immune and cancer biology. In this sense, the IPI dataset and the accompanying curation by these archetypal assignments can now serve as a rich resource to gain deeper understanding of cancer immunity in patients and therefore serve as a framework to direct immunotherapies to the most relevant biology.

Limitation of the Study:

There are a few caveats to note that impose limitations on this study that endeavors to define the dominant tumor immune archetypes across many tumor types. First, we used a single standardized processing approach regardless of the tumor type which was not ideal but necessary due to the limited quantity of available tumor tissue. Nevertheless, both simple immunofluorescence and identification of archetype specific chemokine expression program would appear to support that we reliably recovered cells across tissues, but some archetypes such as the myeloid centric appear to be more impacted (Figure 2H-I) and further investigation of the spatial landscape of the tumor immune archetypes would be needed. In the future, a combination of multiplexed imaging technologies such ion beam imaging (MIBI) [START_REF] Angelo | Multiplexed ion beam imaging of human breast tumors[END_REF][START_REF] Keren | A Structured Tumor-Immune Microenvironment in Triple Negative Breast Cancer Revealed by Multiplexed Ion Beam Imaging[END_REF] and single-cell spatial transcriptomics [START_REF] Asp | A Spatiotemporal Organ-Wide Gene Expression and Cell Atlas of the Developing Human Heart[END_REF][START_REF] Hu | ZipSeq: barcoding for real-time mapping of single cell transcriptomes[END_REF] will help mitigate this issue and enable analysis. In this way, archetype-discovery will play an important role in selecting tissues containing similar biology for studies.

Second, we focused our study on 10 major tumor associated immune cell types across a wide variety of human cancer types, and this potentially biases analyses to dominant tumor archetypes because we only examine pre-defined cell populations. Interestingly, we demonstrated that the archetypes defined by these 10 features often correlate with enrichment of immune cell types that were not part of the clustering (Figure 6A). However, some other populations previously described as important for the tumor immune response, such as neutrophils, γδ T cells [START_REF] Coffelt | IL-17-producing γδ T cells and neutrophils conspire to promote breast cancer metastasis[END_REF][START_REF] Wellenstein | Loss of p53 triggers WNT-dependent systemic inflammation to drive breast cancer metastasis[END_REF] or specific transcriptomic states of cDCs [START_REF] Barry | A natural killer-dendritic cell axis defines checkpoint therapyresponsive tumor microenvironments[END_REF][START_REF] Cheng | A pan-cancer single-cell transcriptional atlas of tumor infiltrating myeloid cells[END_REF][START_REF] Lopes | Distinct metabolic programs established in the thymus control effector functions of γδ T cell subsets in tumor microenvironments[END_REF][START_REF] Maier | A conserved dendritic-cell regulatory program limits antitumour immunity[END_REF] were unable to be analyzed or appeared to be variable within archetypes (Figure 6C andS6). Future work using single-cell omics across tumor types -not only transcriptomic based due to the sensitivity of some of those cell types -will be critical. The combination of single-cell profiling in the context of immune archetypes would determine if the 12 tumor archetypes identified in our study can be further subdivided or if the cell subsets identified using these technologies generally align with our current landscape [START_REF] Cheng | A pan-cancer single-cell transcriptional atlas of tumor infiltrating myeloid cells[END_REF][START_REF] Oliveira | Phenotype, specificity and avidity of antitumour CD8+ T cells in melanoma[END_REF][START_REF] Slyper | A single-cell and single-nucleus RNA-Seq toolbox for fresh and frozen human tumors[END_REF][START_REF] Zhang | Single-Cell Analyses Inform Mechanisms of Myeloid-Targeted Therapies in Colon Cancer[END_REF].

Finally, using the IPI cohort we developed archetype specific gene signatures that enabled us to explore the association with disease outcome on a larger scale in the TCGA cohort.

The two cohorts present important differences in term of demographics, tissue histology and patient clinical history. Therefore, an extension of this work on a bigger cohort with the addition of clinical trials designed to demonstrate the benefits of targeting tumor archetypes would be incredibly important.

STAR METHODS

RESOURCE AVAILABILITY

Lead Contact: Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Matthew F. Krummel (matthew.krummel@ucsf.edu).

Materials availability: This study did not generate new unique reagents.

Data and Code availability:

Single-cell RNA-seq and bulk RNAseq data have been deposited at GEO and are publicly available as of the date of publication. Accession numbers are listed in the key resource table. Microscopy and flow data reported in this paper will be shared by the lead contact upon request.

All original code has been deposited at GitHub and is publicly available as of the date of publication. DOIs are listed in the key resources table.

Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

EXPERIMENTAL MODEL AND SUBJECT DETAILS

Human tumor collection of the UCSF Immunoprofiler Initiative (IPI)-Tumor samples for the Immunoprofiler was transported from various cancer operating rooms (ORs) as well as from outpatient clinics. All patients consented by the UCSF IPI clinical coordinator group for tissue collection under a UCSF IRB approved protocol (UCSF IRB# 20-31740). Samples were obtained after surgical excision with biopsies taken by Pathology Assistants to confirm the presence of tumor cells. Patients were selected without regard to prior treatment. Freshly resected samples were placed in ice-cold PBS or Leibovitz's L-15 medium in a 50 mL conical tube and immediately transported to the laboratory for sample labeling and prepare either the whole tissue for digestion into single-cell suspension or a part of the tissue was sliced and preserved for imaging analysis. All Clinical information on the different patient of the cohort can be found in supplementary table S1.

Assembling Cohorts

UCSF Immunoprofiler initiative (IPI):

An initial set of 427 bulk RNAseq patient samples in the live compartment were evaluated based on an in-house metric, the EHK score, that serves as a measure of data quality. Each sample is given a score of 0 through 10 depending on the number of EHK genes that are expressed above a precalculated minimum threshold. The threshold was learned from our data by examining the expression distributions of EHK genes and validated using the corresponding distributions in TCGA. A score of 10 represents the highest quality data where 10 out of 10 EHK genes are expressed above the minimum threshold. Filtering for samples with an EHK score of EHK8, EHK9 and EHK10 reduced the sample set to 298 patient samples. The sample set was then filtered to remove all adjacent normal samples and all biological replicates, reducing the sample set further to 260 samples. These 260 patient samples are the IPI cohort along with 199 overlapping tcell compartment samples and 189 myeloid compartment samples.

TCGA: Tumor RNAseq counts and TPM along with curated clinical data for 13 cancer types (Bladder urothelial carcinoma (BLCA), Colon adenocarcinoma (COAD), glioblastoma multiforme (GBM), Head and Neck squamous cell carcinoma (HNSC), Kidney renal clear cell carcinoma (KIRC), Liver hepatocellular carcinoma (LIHC), Lung adenocarcinoma (LUAD), Ovarian serous cystadenocarcinoma (OV), Pancreatic adenocarcinoma (PAAD), Sarcoma (SARC), Skin Cutaneous Melanoma (SKCM), Uterine Corpus Endometrial Carcinoma (UCEC), Uterine Carcinosarcoma (UCS)), from the Toil recompute [START_REF] Vivian | Toil enables reproducible, open source, big biomedical data analyses[END_REF] data in the TCGA Pan-Cancer (PANCAN) cohort, were downloaded from the UCSC Xena browser [START_REF] Goldman | Visualizing and interpreting cancer genomics data via the Xena platform[END_REF]. The initial set of 4677 tumor samples was filtered down to include primary solid tumors and metastatic sample types (sample type codes = 01 & 06) only, to parallel the IPI cohort sample types as accurately as possible. This reduced the patient sample set to 4341 tumor samples.

METHOD DETAILS

Human tissue digestion and Multiparametric Flow Cytometry staining and sorting-Tumor or metastatic tissue was thoroughly chopped with surgical scissors and transferred to GentleMACs C Tubes (Miltenyi Biotec) containing 20 uL/mL Liberase TL (5 mg/ml, Roche) and 50 U/ml DNAse I (Roche) in RPMI 1640 per 0.3 g tissue. GentleMACs C Tubes were then installed onto the GentleMACs Octo Dissociator (Miltenyi Biotec) and incubated for 45min according to the manufacturer's instructions. Samples were then quenched with 15 mL of sort buffer (PBS/2% FCS/2mM EDTA), filtered through 100 um filters and spun down. Red blood cell lysis was performed with 175 mM ammonium chloride if needed.

Cells were then incubated with Human FcX (Biolegend) to prevent non-specific antibody binding. Cells were then washed in DPBS and incubated with Zombie Aqua Fixable Viability Dye (Thermo). Following viability dye, cells were washed with sort buffer and incubated with cell surface antibodies mix diluted in the BV stain buffer (BD Biosciences) following manufacturer instruction for 30 minutes on ice in the dark and subsequently fixed in either Fixation Buffer (BD Biosciences) or in Foxp3/Transcription Factor Staining Buffer Set (eBioscience) if intracellular staining was required. For cell sorting single cell suspension was stained with cell surface antibodies mix diluted in BV stain buffer (BD Biosciences) following manufacturer instruction for 30 minutes on ice in the dark and subsequently washed three time and resuspended in of sort buffer (PBS/2% FCS/2mM EDTA) after filtering in 100um filter (thermo). When possible, a maximum of 6 different tumor associated cell subsets were sorted with a BD FACSAria Fusion into 1.5mL Eppendorf tubes containing 150ul of lysis buffer used for mRNA extraction (Invitrogen). The 6 different tumor associated cell types between 5000 and 50000 cells were sorted based their protein expression profile as describe here: 'live cells' all cells negative for the viability dye; Tconv: negative for viability dye CD45+ CD3+ CD19-CD56-CD25-; Treg: negative for viability dye CD45+ CD3+ CD19-CD56-CD4+ CD25+ ; Myeloid: negative for viability dye CD45+ CD3-CD19-CD56-HLADR+; CD90+ CD44+ Stroma: negative for viability dye CD45-CD44+ CD90+; Tumor: negative for viability dye CD45-CD44-CD90-. Names of the different sorted populations were given based on the presumptive populations, and we accepted that both due to slight sorting contamination and due to the limits of antibodies to completely 'define' a population, that other populations would occasionally infiltrate these e.g., some 'activated' CD4 cells express CD25 and might infiltrate this RNA sample.

Cell sorting for bulk RNA sequencing

Cell sorting for single-cell RNA sequencing of the tumor associated mononuclear phagocytic cells.-For the discovery cohort scRNA-seq, live CD3-CD19/20-CD56-SSC-A dim CD16 dim (to exclude neutrophils) cells were sorted from a melanoma involved draining LN on a BD FACSAria Fusion. For the discovery cohort the same sorting strategy was applied to melanoma, kidney and head and Neck primary tumors. After sorting, cells were pelleted and resuspended at 1x10e3 cells/ul in 0.04%BSA/PBA and loaded onto the Chromium Controller (10X Genomics). Samples were processed for single-cell encapsulation and cDNA library generation using the Chromium Single Cell 3' v2 Reagent Kits (10X Genomics). The library was subsequently sequenced on an Illumina HiSeq 4000 (Illumina). All samples were sequenced at 25,000 reads per cell. Filters used for specific fluorophores are: Spectral Gold (Semrock) filter was used for Rhodamine 6G and 87 HE DAPI (Zeiss) for DAPI. H&E slides were scanned in brightfield with the same objective and a HV-F202SCL CCD camera (Hitachi).

Bulk RNA-Sequencing

Library Preparation: The mRNA was isolated via DynaBeads Direct, then converted into amplified cDNA using Tecan Ovation RNA-Seq System V2 kit, following the manufacturer protocol. The dsDNA went through tagmentation, amplification and clean up with AMPure XP beads steps, using the Illumina Nextera XT DNA Library Prep Kit. Quality control analysis was performed on the resulting pooled libraries with the Agilent Bioanalyzer HS DNA chip to assess fragment size distribution and concentration.

Sequencing:

The pooled libraries were then sequenced via single-read MiSeqMiniSeq to check if more than 10 percent of the reads aligned to coding regions and contained more than 1,000 unique reads in total. Sequencing efficiency was determined by protein-coding read fractions. The protein-coding read fractions were then used to normalize libraries for high depth sequencing with the goal of obtaining an average of 10M protein-coding reads per sample. These libraries were submitted the UCSF Center for Advanced Technology for 100bp paired end (PE100) sequencing on the HiSeq4000.

Bioinformatic Data Processing

Bulk RNAseq: The RNAseq reads were first aligned to reference of rRNA and mitochondrial rRNA using BWA-mem to deplete the dataset of any remaining rRNA. The remaining reads were aligned to the Ensembl GRCh38.85 transcriptome build using STAR [START_REF] Dobin | STAR: ultrafast universal RNA-seq aligner[END_REF]. Gene expression was computed from the alignments in counts and TPM (transcripts per million) using RSEM [START_REF] Li | RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome[END_REF].

QC of Bulk RNAseq: QC of the of the RNAseq was done to identify any swapped samples. DEG was performed on all pairs of the tcell, myeloid, treg, stroma and tumor compartments. The top 500 genes by p value were identified for each DEG pair giving 5000 genes which were reduced to 2171 genes after removing redundant genes. Principal component analysis (PCA) was done on the mean-centered and scaled logCPM (log counts per million) of the counts of these genes. K-means clustering was performed on the first 4 principal components with K=5. This yielded 5 clusters that were mostly made up of one compartment. Any samples that clustered with a different compartment were considers swaps except for the tumor compartment which we assumed would be promiscuous.

Data pre-processing of 10x Genomics Chromium scRNA-seq data: Feature-barcode matrices were obtained for sample by aligning the raw fastqs to GRCh38 reference genome (annotated with Ensembl v85) using the Cellranger count. Raw feature-barcode matrices were loaded into Seurat 3.1.5 [START_REF] Stuart | Comprehensive Integration of Single-Cell Data[END_REF] and genes with fewer than 3 UMIs were dropped from the analyses. Matrices were further filtered to remove events with greater than 20% percent mitochondrial content, events with greater than 50% ribosomal content, or events with fewer than 100 total genes. The cell cycle state of each cell was assessed using a published set of genes associated with various stages of human mitosis [START_REF] Dominguez | A high-resolution transcriptome map of cell cycle reveals novel connections between periodic genes and cancer[END_REF].

Data quality control and Normalization for scRNA-seq data:

The filtered count matrices were normalized, and variance stabilized using negative binomial regression via the scTransform method offered by Seurat [START_REF] Hafemeister | Normalization and variance stabilization of single-cell RNA-seq data using regularized negative binomial regression[END_REF]. The effects of mitochondrial content, ribosomal content, and cell cycle state were regressed out of the normalized data to prevent any confounding signal. The normalized matrices were reduced to a lower dimension using Principal Component Analyses (PCA) and the first 30 principal coordinates per sample were subjected to a non-linear dimensionality reduction using UMAP reduction. Clusters of cells sharing similar transcriptomic signal were identified using the Louvain algorithm, and clustering resolutions varied between 0. 6 and 1.2 based on the number and variety of cells obtained in the datasets.

Data integration and Batch correction for the scRNA seq validation dataset:

The individual processed objects per library were normalized, and variance stabilized using negative binomial regression via the scTransform method offered by Seurat [START_REF] Hafemeister | Normalization and variance stabilization of single-cell RNA-seq data using regularized negative binomial regression[END_REF]. Counts matrices were merged into a single Seurat object and the batch (or library) of origin was stored in the metadata of the object. The log-normalized counts were reduced to a lower dimension using PCA and the individual libraries were aligned in the shared PCA space in a batch-aware manner (Each individual library was considered a batch) using the Harmony algorithm [START_REF] Korsunsky | Fast, sensitive and accurate integration of single-cell data with Harmony[END_REF]. The resulting Harmony components were used to generate a batch corrected UMAP, and to identify clusters of transcriptionally similar cells across each of the 11 samples. The 5 Genes signature identified in the discovery dataset from the melanoma tumor sample was used to identify monocytes, macrophages, cDC2 and cDC1 respectively.

Single Cell RNAseq intra-sample heterotypic doublet detection: All libraries were further processed to identify heterotypic doublets arising from the 10X sample loading. Processed, annotated Seurat objects were processed using the DoubletFinder package [START_REF] Mcginnis | DoubletFinder: Doublet Detection in Single-Cell RNA Sequencing Data Using Artificial Nearest Neighbors[END_REF]. Briefly, the cells from the object are modified to generate artificial duplicates, and true doublets in the dataset are identified based on similarity to the artificial doublets in the modified gene space. The prior doublet rate per library was approximated using the information provided in the 10x knowledgebase (https://kb.10xgenomics.com/hc/ en-us/articles/360001378811) and this was corrected to account for homotypic doublets using the per-cluster numbers in each dataset.

Imaging bioinformatics: Each image was processed and analyzed in Imaris v9.5 (Bitplane Inc.). We first applied a median filter with size 3x3x1 to every channel of the image. The background of each channel was then subtracted using a filter width of 100um. Next, we used the Spots function to detect cell nuclei in the DAPI channel. We then applied a binary threshold on the max intensity of the CD45 channel to classify each spot as positive or negative for CD45.

We then cropped out a 650um x 650um section containing the tumor border from each histological image. This same section of tissue was found in the immunofluorescence image by looking for tissue features shared between the Immunofluorescence (IF) and histological images. Using the CD45 counts described above, two representative images from each immune class were selected.

QUANTIFICATION & STATISTICAL ANALYSIS

Analysis using 3 Features -Tcell, Myeloid, CD90+ CD44+ Stroma

Gene Signature: Gene signatures for each of the features (Supplementary Table 4) were generated by performing differential gene expression (DEG) analysis, between the different sorted RNAseq compartments, using edgeR (Robinson et al., 2010), limma [START_REF] Ritchie | limma powers differential expression analyses for RNA-sequencing and microarray studies[END_REF] and Voom [START_REF] Law | voom: precision weights unlock linear model analysis tools for RNA-seq read counts[END_REF], a function of limma that modifies RNAseq data for use with limma. Expression counts for all EHK10 samples in each of the compartments (tcell, myeloid, stroma, treg and tumor) were selected to ensure the DEG analysis was done using the highest quality data. The median adjusted p-value for these DEGs were values much less than 1.00 x 10 -55 across all comparisons. For the Tcell feature gene signature, DEG analysis was done between the tcell compartment and the stroma, myeloid, and tumor compartments separately. The intersection of the top 50 genes by log fold change (logFC) for each of these three comparisons was designated the Tcell feature gene signature. Similarly, a CD90+ CD44+ Stroma feature gene signature was designated as the intersection between DEG comparisons of the stroma compartment versus the myeloid, tumor and tcell compartments.

For the Myeloid feature gene signature, it was necessary to use the top 100 genes by logFC from the DEG between the myeloid compartment and tcell compartment to compensate for the fact that myeloid cells are more auto fluorescent and thus more likely to contaminate the tcell compartment sort. The stroma and tumor compartments are sorted from CD45-and are less likely to experience this issue. The process of feature gene signature selection for each of the three features was visualized with volcano plots that showed the DEG results and Venn diagrams that showed the intersection of the top genes from each DEG and the resulting feature gene signature.

Flow Score: The flow scores for the flow cytometry population ratios were calculated by scaling the ratios from 1-100. This was done by calculating the percentile rank of each population ratio relative to the other population ratios. For example, a percentile rank of 65 means that that 65% of the population ratios are below that ratio. This array of percentile ranks was then further transformed to percentile ranks to increase the spread in the scores. Boxplots were generated to visualize the scores by cancer type using Seaborn a python plotting module.

Gene Signature Score: The feature gene signature scores were calculated using an m x n matrix where m represented the TMM [START_REF] Robinson | A scaling normalization method for differential expression analysis of RNA-seq data[END_REF] normalized logCPM (log2 counts per million) expression of the feature signature genes and n represented the selected sample set. TMM normalization is a method that adjusts library size so samples can be compared, since sequencing depth can differ between samples.

The expression of each gene was converted to percentile ranks across the samples using the SciPy [START_REF] Virtanen | SciPy 1.0: fundamental algorithms for scientific computing in Python[END_REF] Python module.

E P

where: E = m x n matrix of gene expression P = m x n matrix of gene expression percentile ranks across samples A score was generated for each sample n follows:

Score n = percentile(Average(P n ))

where: P n = column n of P corresponding to sample n (The second percentile transformation was employed to increase the spread between the scores). Boxplots were generated to visualize the distribution of feature scores by cancer type with statistical analysis between archetypes performed using the statannot python package. We used the Mann-Whitney test between two archetypes, as we could not assume normal distributions, with the Bonferroni correction for multiple comparisons. The significance levels were annotated on the boxplots using the P-Value annotation legend below. For the IPI gene signature scores for the Tcell, Myeloid, CD90+ CD44+ Stroma were calculated using the live compartment. For the TCGA cohort the gene signature scores for the Tcell, Myeloid and CD90+ CD44+ Stroma features were calculated for all tumor type.

Cross-Whisker plots, that show the median value, the interquartile range (IQR) in both the x and y directions and the identity line, were made to compare the median Tcell, Myeloid and CD90+ CD44+ Stroma feature scores in the IPI cohort with the TCGA cohort. Based on the Cross-Whisker plots only the BLCA, COAD, UCS, UCEC, OV, HNSC, KIRC, SARC and SKCM, cancer type were used for further analysis as their, median scores per tumor type appeared to correlate well with the IPI patient samples. The rationale behind including only tumor type that correlated well with IPI was to ensure that further analysis in TCGA would have extrapolated relevance in the IPI patient samples. The excluded cancer type lack of correlation with IPI patient samples could be attributed to variety of reasons that cannot be easily controlled, e.g. technical variability in upstream processing steps and cohort race differences due to local demographic biases.

Gene Signature Score Validation: Based on the assumption that the flow cytometry population ratios represent the true feature abundance in the patient sample, correlation between flow cytometry population ratios and feature gene signatures were used to validate the feature gene signatures. On average 9 to 10 samples per cancer type, that had both flow cytometry data and RNAseq data in the relevant compartments, were selected (Supplementary Table 3). There were fewer samples selected in PDAC (Pancreatic), PNET(Neuroendocrine) and GBM(Glioblastoma) cancer type as there were not many representative samples for theses cancer in the IPI cohort. Correlation plots were generated using the scores and their corresponding flow population ratios and the Spearman's rho was calculated for each of the correlations.

Additionally, the Tcell and Myeloid feature samples were submitted to CIBERSORT [START_REF] Newman | Robust enumeration of cell subsets from tissue expression profiles[END_REF]. The TPMs of all protein coding genes of the samples were input as the mixture file, the built in LM22 gene signatures were used and quantile normalization was disabled as this was the recommended setting for RNAseq data. The CIBERSORT score used for the Tcell feature was the summation of the T cells CD8, T cells CD4 naïve, T cells CD4 memory resting, T cells CD4 memory activated, T cells follicular helper, T cells regulatory (Tregs) and T cells gamma delta cell type relative fractions. The CIBERSORT score used for the Myeloid feature was the summation of the Macrophages M0, Macrophages M1, Macrophages M2, Dendritic cells resting, Dendritic cells activated, Mast cells resting, Mast cells activated, Eosinophils and Neutrophils cell type relative fractions. Correlation plots were generated using these CIBERSORT scores and their corresponding flow population ratios and the Spearman's rho was calculated for each of the correlations.

Clustering:

The feature gene signature scores in the IPI cohort were clustered using SCANPY [START_REF] Wolf | SCANPY : large-scale single-cell gene expression data analysis[END_REF] a Python-based toolkit. The feature clustering was done in three steps. In the first step, a neighborhood graph was constructed using the k-nearest neighbors (KNN) algorithm with Euclidean distances. In the second step, the neighborhood graph was clustered using the Louvain method [START_REF] Blondel | Fast unfolding of communities in large networks[END_REF], a community detection algorithm that maximizes network modularity. In the third step, the resulting clustering was evaluated using the Davies Bouldin Index (DBI) (Davies and Bouldin, 1979), a metric that assesses the ratio of the intra-cluster distance to inter-cluster distance. The lower the DBI the better the separation between clusters, the more compact the samples within clusters and hence the better the clustering. The clustering steps were repeated incrementally over values of n between 3 to 300 for the nearest neighbors for KNN. At each iteration of n, the values for the resolution parameter for Louvain clustering were varied between 0.3 to 2 in 0.1 increments. The formula below was used to optimize the nearest neighbor and resolution parameters for final cluster selection. Where: n = [3,4,5…,300] (KNN nearest neighbor) res = [03,0.4,0.5…,2] (resolution)

In the TCGA cohort the feature signature score clustering was done as above except for the values for the range of n for KNN which was 3-1,500 as the number of samples in the TCGA cohort were greater.

The final clustering selections were given archetype labels based on the prevalence and distributions of the features within the cluster which were visualized as violin plots of all clusters in each feature. Both the IPI and TCGA final clustering selections had six clusters and hence six archetypes. The final clustering selection (IPI-was visualized as UMAP [START_REF] Mcinnes | UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction[END_REF], a dimensionality reduced projection of the clustering. The min-dist UMAP parameter, which is the minimum distance between points, controls how close together points are placed in the low-dimensional space of the UMAP was set as of 0.25. Population flow proportions and phenotype gene signatures were overlayed on the UMAP clustering.

Univariate survival analysis on the TCGA data was performed per cancer type and multivariate survival regression was performed across all cancer types using cancer types as a covariates using lifelines. Pie charts were made to compare the archetype abundance per tumor type between IPI and TCGA. Additionally, 3D plots of IPI score and TCGA score visualized the distribution of the unclustered scores.

Phenotype Gene Signatures:

The Chemokine gene signature (Supplementary Table 5) was a selection of 39 chemokines [START_REF] Nagarsheth | Chemokines in the cancer microenvironment and their relevance in cancer immunotherapy[END_REF] of interest and was evaluated in the live compartment to generate a score. Hierarchical clustered heatmaps were made to assess chemokine expression in archetypes. The median TPM expression per archetype of each of the 39 chemokines in the Chemokine gene signature were clustered using the Euclidean distance metric and the Ward linkage method. A bubble plot was made to compare the median TPM expression of Chemokine gene signature genes in IPI and TCGA in each archetype. The bubble size represented values of the median TPM in each archetype transformed to a z-score.

Analysis using 6 Features -Tcell, Myeloid, CD90+ CD44+ Stroma, CD4, CD8, Tregs Gene Signature -CD4, CD8 and Tregs: CD4 high and CD8 high samples were identified based on the flow proportion of (CD4+, CD25-, FOXP3-) / CD3+ and (CD8+, CD4-) / CD3+ respectively. There were nine samples that were both EHK10 in the tcell compartment of RNAseq and had a (CD4+, CD25-, FOXP3-) / CD3+ flow proportion greater than 60%, these were labeled CD4 high. Similarly, there were 12 samples that were EHK10 in the tcell compartment of RNAseq that had a (CD8+, CD4-) / CD3+ flow proportion greater than 60%, these were labeled CD8 high. The CD4 gene signatures were generated by performing DEG analyses using edgeR and voom on the tcell compartment counts of CD4 high samples versus CD8 high samples. The top 50 genes by adjusted P-value were selected and all genes within this gene set that had a positive logFC were designated the CD4 feature gene signature. The same was done for the CD8 feature gene signature except the DEG was CD8 high samples versus CD4 high samples (Supplementary Table 4). A hierarchical clustered heatmap was made using the CD4 and CD8 gene signature logCPM values for all tcell compartment patient samples. The logCPM values were mean-centered and scaled prior to clustering using the Euclidean distance metric and the Ward linkage method.

The Treg feature gene signature followed all the same steps as the Tcell feature gene signature but with four DEG comparisons instead of three, the treg compartment versus the tcell, myeloid, stroma and tumor compartments. The resulting intersection of these four comparisons was designated the Treg feature gene signature.

Gene Signature Score: Gene signature scores for the three new features were calculated as described above for three features. The CD4 and CD8 gene signatures were evaluated in the tcell compartment and the Treg gene signature was evaluated in the live compartment. Correlation between flow cytometry population ratios and feature gene signatures were plotted used to validate the feature gene signatures.

Calculated Score for Missing Data: The IPI cohort consists of 260 patient samples that all have a live compartment sequenced. However, there are only 199 overlapping patient samples with the tcell compartment sequenced. Since the CD4 and CS8 gene signature score were evaluated in the tcell compartment a need for a method to calculate scores for missing data was warranted. Gene signature scores were calculated for missing patient samples by leveraging the fact that there was very high correlation between gene signatures scores and their corresponding flow population ratios. For example, the correlation between the Tcell gene signature score and the CD3+ / live population ratio from flow cytometry had a Spearman's rho of 0.91. The missing scores were calculated by modeling the correlation using linear regression and using the regression equation and the flow population ratios to calculate the score. Of the 61 missing tcell compartment samples 5 did not have flow data.

Clusterin:

The feature gene signature scores for the six features were clustered as described above. The final clustering selections had eight clusters and hence eight archetypes. The final clustering was visualized as a UMAP. Population flow proportions and phenotype gene signatures were overlayed on the UMAP clustering.

Phenotype Gene Signatures: Hierarchical clustered heatmaps were made to assess chemokine expression in archetypes. The median TPM expression per archetype of each of the 39 chemokines in the Chemokine gene signature were clustered using the Euclidean distance metric and the Ward linkage method.

Analysis using 10 Features -Tcell, Myeloid, CD90+ CD44+ Stroma, CD4, CD8, Tregs, Macrophages, Monocytes, cDC1, cDC2

Gene Signature -Macrophages, Monocytes, cDC1 and cDC2: Gene signatures were generated from the discovery single cell RNAseq on tumor associated mononuclear phagocytes (MNP) dataset using a melanoma tumor sample. Gene signatures were generated using the FindMarkers/FindAllMarkers functions in Seurat. The genes that were considered significant had a log-fold change greater than 0.4, an adjusted p value less than or equal to 0.05 (based on Bonferroni correction) and were expressed in at least 35% of the groups.

Cluster marker genes were identified as the top upregulated genes of the DEGs between clusters and curated based on their log-fold change as well as low expression by the other clusters. The different subsets of mononuclear phagocytic cells were identified by comparing cluster marker genes with public sources referenced in the text. The list of the different genes for each subset is listed in (Supplementary Table 4)

Gene Signature Score: Gene signature scores for the four new features were calculated as described above for three features. The gene signatures were evaluated in the myeloid compartment Correlation between flow cytometry population ratios and feature gene signatures were plotted used to validate the feature gene signatures.

Calculated Score for Missing Data: Of the of 260 patient samples that all have a live compartment sequenced there were only 189 overlapping patient samples with the myeloid compartment sequenced. The same method used to calculate missing data scores for the CD4 and CD8 gene signatures was used for missing data in the Macrophage, Monocyte, cDC1 and cDC2 gene signature. Of the 71 missing myeloid compartment samples 15 did not have flow data.

Clustering:

The feature gene signature scores for the 10 features were clustered as described above. The final clustering selections had 12 clusters and hence 12 archetypes. The final clustering was visualized as a UMAP. Population flow proportions and phenotype gene signatures were overlayed on the UMAP clustering. Alluvial plots were made using RAWgraphs [START_REF] Mauri | RAWGraphs: A Visualisation Platform to Create Open Outputs[END_REF] to visualize the stability of clusters as the clustering progressed from six features to ten features.

Phenotype Gene Signatures: Phenotype gene signatures (Supplementary Table 5) were used to further characterize the immune archetypes obtained by clustering feature scores generated from the feature gene signatures. The phenotype gene signatures were assembled from standard gene signatures, derived from analysis of IPI samples or curated from gene sets found in literature.

The Tcell Exhaustion gene signature was generated by correlating the expression of all protein coding genes in EHK10 samples from the tcell compartment with the expression of five genes, CTLA4, PDCD1, CD38, HAVCR2 and LAG3, a subset of a published T cell exhaustion gene signature (B et al., 2018). The top 50 genes by Spearman's rho for each of the five correlations were selected as gene sets and genes that fell in the intersection of at least four of these gene sets were designated as the Tcell Exhaustion feature gene signature. The Spearman's rho of the CTLA4, PDCD1, CD38, HAVCR2 and LAG3 against the Tcell Exhaustion gene signature genes with additional known exhaustion genes were plotted on heatmap to illustrate the gene signature selection criteria.

The ISG, Senescence, EMT, Fibrosis, Cell Stress DNA damage, Cell Cycle G1-S and Cell Cycle G2-M genes signatures were from literature (Kinker et al., 2019, munoz dp et al, wiley cd et al). A bubble plot was made to visualize the median TPM expression per archetype in the tumor compartment. The bubble size and color represented values of the median TPM in each archetype, transformed to a z-score.

The NK cell gene signature was from literature [START_REF] Barry | A natural killer-dendritic cell axis defines checkpoint therapyresponsive tumor microenvironments[END_REF], the B cells gene signature was from [START_REF] Chen | Single-cell transcriptome and antigen-immunoglobin analysis reveals the diversity of B cells in non-small cell lung cancer[END_REF][START_REF] Van Galen | Single-Cell RNA-Seq Reveals AML Hierarchies Relevant to Disease Progression and Immunity[END_REF] , the Plasma cells gene signature was from [START_REF] Chen | Single-cell transcriptome and antigen-immunoglobin analysis reveals the diversity of B cells in non-small cell lung cancer[END_REF][START_REF] Van Galen | Single-Cell RNA-Seq Reveals AML Hierarchies Relevant to Disease Progression and Immunity[END_REF] and the Mast cell signature was from [START_REF] Cheng | A pan-cancer single-cell transcriptional atlas of tumor infiltrating myeloid cells[END_REF]. A bubble plot was made to visualize the median TPM expression per archetype in the live compartment. The bubble size and color represented values of the median TPM in each archetype, transformed to a z-score.

The Th1, Th2, Th17, Trm and Tex gene signatures were from [START_REF] Kumar | Human Tissue-Resident Memory T Cells Are Defined by Core Transcriptional and Functional Signatures in Lymphoid and Mucosal Sites[END_REF][START_REF] Savas | Single-cell profiling of breast cancer T cells reveals a tissueresident memory subset associated with improved prognosis[END_REF][START_REF] Zhou | Plasticity of CD4+ T cell lineage differentiation[END_REF] [START_REF] Bengsch | Epigenomic-Guided Mass Cytometry Profiling Reveals Disease-Specific Features of Exhausted CD8 T Cells[END_REF].

The Resting Tregs, Suppressive Tregs, Tissue Homing Tregs, Treg Activation, Treg Cytokines, Treg Metabolism were from (Arce [START_REF] Vargas | Fc Effector Function Contributes to the Activity of Human Anti-CTLA-4 Antibodies[END_REF][START_REF] Ephrem | Modulation of Treg cells/T effector function by GITR signaling is context-dependent[END_REF][START_REF] Plitas | Regulatory T Cells Exhibit Distinct Features in Human Breast Cancer[END_REF][START_REF] Zemmour | Single-cell gene expression reveals a landscape of regulatory T cell phenotypes shaped by the TCR[END_REF] The M1, M2 and Costimulatory molecule gene signatures were from [START_REF] Biswas | Tumor-associated macrophages: functional diversity, clinical significance, and open questions[END_REF][START_REF] Cassetta | Human Tumor-Associated Macrophage and Monocyte Transcriptional Landscapes Reveal Cancer-Specific Reprogramming, Biomarkers, and Therapeutic Targets[END_REF][START_REF] Maier | A conserved dendritic-cell regulatory program limits antitumour immunity[END_REF][START_REF] Roberts | Critical Role for CD103(+)/CD141(+) Dendritic Cells Bearing CCR7 for Tumor Antigen Trafficking and Priming of T Cell Immunity in Melanoma[END_REF]. A bubble plot was made to visualize the median TPM expression per archetype in the myeloid. The bubble size and color represented values of the median TPM in each archetype, transformed to a z-score.

The MHC Class I gene signature was from HGNC ( https://www.genenames.org). The ISG gene signature was from [START_REF] Combes | Global absence and targeting of protective immune states in severe COVID-19[END_REF]. A box plot was made to visualize the scores per archetype that were calculated in the live compartment.

A hierarchical clustered heatmap was made to visualize the Th1, Th2, Th17, Trm, Tex, Resting Tregs, Suppressive Tregs, Tissue Homing Tregs, Treg Cytokines, Treg Activation, Treg Metabolism, NK cells, Plasma cells, B cells, M1, M2, Costimulatory molecule, and ISG gene signature gene expression. These genes for each of the phenotype gene signatures were assessed different RNAseq compartments as described above. The median TPM expression per archetype of each of the genes in the gene signatures were clustered using correlation distance metric and the average linkage method.

Immune Archetype Gene Signatures

Live: Archetype gene signatures were generated by DEG analysis, between the live compartment counts of samples in each of the 12 archetypes, using limma and Voom. The intersection between the top 3000 genes by logFC of each of 11 DEGs per archetype was assigned as an initial gene signature. If the initial gene signatures had less than 20 genes it was designated as the archetype gene signature. Otherwise, the archetype gene signature was designated as the top 20 genes with both the lowest coefficient of variation (CV) of the log10(TPM+0.001) expression and non-zero expression in at least 80% of the samples in the archetype. A hierarchical clustered heatmap was made to visualize archetype gene signature gene expression. The median TPM expression per archetype of each of the genes in the gene signatures were clustered using correlation distance metric and the average linkage method.

Tumor: Archetype gene signatures were generated by DEG analysis, between the combined tumor and epcam compartment counts of samples in each of the 12 archetypes, using edgeR and Voom. The intersection between the top 100 genes by logFC of each of 11 DEGs per archetype was assigned as an initial gene signature. If the initial gene signatures had less than 10 genes it was designated as the archetype gene signature. Otherwise, the archetype gene signature was designated as the top 10 genes with both the lowest coefficient of variation (CV) of the log10(TPM+0.001) expression and non-zero expression in at least 80% of the samples in the archetype. A hierarchical clustered heatmap using correlation distance metric and the average linkage method was made for the median TPM expression per archetype of each of the genes in the gene signatures.

Live Archetype Gene Signature Score: Using TCGA data, a gene signature score for each archetype was calculated using the gene signature score method above. Each sample has 12 scores for each archetype. Each sample was assigned to the archetype for which it had the highest score and if the highest score was tied between archetypes, the sample was excluded from the analysis.

Univariate survival analysis was performed per cancer type and multivariate survival regression was performed across all cancer types using cancer types as a covariates. The analysis was done using lifelines (Davidson-Pilon), a survival analysis Python library. If an archetype had less than fifteen representative samples, it was excluded from analysis. Pie charts were made to compare the archetype abundance per cancer type between IPI and TCGA.

ADDITIONAL RESSOURCES

To facilitate usage of our data for the wide research community, we will provide access to transcriptomic and compositional data for each archetype: https://datalibrary.ucsf.edu/ public-resources.

Supplementary Material

Refer to Web version on PubMed Central for supplementary material.

• UCSF Immunoprofiler Initiative is a set of human tumors with multimodal linked data.

•

Clustering upon 10 features identifies 12 unique tumor Archetypes spanning cancer types.

•

Each Archetype concentrates similarities in additional immune and tumor features.

• Dominant archetypes aid in tumor classification and identifying therapeutic targets.

Tumors can be extremely heterogenous across multiple parameters. By identifying and classifying recurrent immune features across 12 different tumor types, a framework is provided towards understanding tumor immunity as well as identifying broad as well as common therapeutic targets.

  min[DBI(Cluster n, res )]

Figure 1 :

 1 Figure 1: Generation and validation steps of T cells, Myeloid cells and Stromal cells features from solid tumors using flow cytometry and bulk RNA-sequencing. A. Details of the Immunoprofiler initiative (IPI) cohort tumor samples collection, colorcoded by anatomical region and annotated with case numbers of bulk-RNA sequencing of viable cells sorted from fresh surgical tumor specimens and total number of samples with flow cytometry data. B. Description of the processing pipeline for digesting fresh tumor specimens into single cell suspension, submitting to multi-parametric flow cytometry for immune phenotyping and cell sorting into six different cell population compartments (live (viable cells), tcell (conventional T cells), treg (T regulatory cells), myeloid (myeloid cells, stroma (CD90+ CD44+ stromal cells), and tumor (tumor cells). C. Box and whisker plots of flow score for T cells (n=200), Mononuclear phagocytes (n=159), and Stromal cells (n=121) based on population percent in tumor specimens measured by flow cytometry (see Supplementary table S2 for details by cancer type). D.-Left-Gene score calculation method

Figure 2 :

 2 Figure 2: Identification of coarse immune archetypes in solid tumors using Louvain clustering on two independent datasets. A,E. UMAP display using KNN and Louvain clustering of tumor immune archetypes using Tcell, Myeloid and CD90+ CD44+ Stroma features to cluster patients in the IPI (A) and TCGA (E) cohorts. Each dot represents a single patient. B,F. UMAP overlays of the Tcell, Myeloid and CD90+ CD44+ Stroma features in the IPI (B) and TCGA (F) cohorts. C. Violin plots of Tcell, Myeloid and CD90+ CD44+ Stroma features for each cluster/archetype in IPI cohort. D. Table summarizing the six cluster/archetypes with descriptions based on the level of the Tcell, Myeloid and CD90+ CD44+ Stroma features. G,H. Representative Immunofluorescence of tumor specimens using CD45 (red) and DAPI

Figure 3 :

 3 Figure 3: Coarse immune archetypes are independent of tissue origin and associated to overall survival. A. Left-UMAP display, and graph-based clustering of 3-feature tumor immune archetypes color-coded by tumor type-Right stacked bar plot of the tumor type distribution for 3-feature archetypes. B. Kaplan-Meier overall survival curves for each immune tumor archetype identified in the TCGA cohort. C. (Left) Pie charts representing distributions of each archetype by cancer type in the IPI (top) and TCGA (bottom) cohorts. (Right) Kaplan-Meier overall survival curve for each immune tumor archetype identified in the TCGA cohort for Kidney renal clear cell carcinoma (KIRC), Skin cutaneous melanoma (SKCM), Sarcoma

Figure 4 :

 4 Figure 4: Inclusion of T cell subset features subdivide immune archetypes by CD4 to CD8 ratio. A. Box and whisker plots of feature gene signature scores for CD4+ regulatory T cells (Treg feature) out of the live compartment, CD4+ and CD8+ conventional T cells (CD4 and CD8 features) out of the tcell compartment of patients in the IPI cohort. B. Heatmap and hierarchical clustering of CD4 (yellow) and CD8 (blue) feature genes' normalized expression, for patients in the tcell compartment. C. (Left) UMAP display and graph-based clustering of tumor immune archetypes using Tcell, Myeloid, CD90+ CD44+ Stroma, CD4, CD8 and Treg features to cluster patients in the IPI cohort. Each dot represents a single patient. (Right) Table summarizing the eight cluster/archetypes with descriptions based on the abundance of the Tcell, Myeloid CD90+ CD44+ Stroma, CD4, CD8 and Treg features. D. Box and whisker plot of log2 CD4+ to CD8+ conventional T cell feature gene signature score ratio in tumor for each of the clusters/archetypes identified with 6-feature clustering. E. Alluvial plot depicting how cluster/archetype membership perpetuates or subdivides from 3 to 6-feature clustering. F. Heatmap and hierarchical clustering of the median chemokine gene expression for each cluster/archetype identified in the 6-feature clustering. G. Box and whisker plot of log2 cDC2 to cDC1 ratio (top) and Mono to Macs ratio (bottom) measured by flow cytometry for each cluster/archetype identified.

Figure 5 :

 5 Figure 5: Single-cell RNA sequencing-derived myeloid signatures refines immune archetypes. A. (Left)Box and whisker plots for the Macrophages, Monocytes, cDC1 and cDC2 features, calculated in the myeloid compartment, from the IPI cohort. (Right) UMAP display and graph-based clustering of tumor immune archetypes using Tcell, Myeloid, CD90+ CD44+ Stroma, Treg, CD4, CD8, Macrophages, Monocytes, cDC1 and cDC2 features to cluster patients in the IPI cohort. B. Alluvial plot depicting how cluster/archetype membership perpetuates or subdivides from 6 to 10-feature clustering. C. (Left) Schematic of a "phylogeny" of the cluster/archetypes as they progressed from 3-feature to 6-feature to 10-feature clustering. D-N. Scatter plots of different features defining the twelve clusters/ archetypes identified in the IPI cohort. O, R. UMAP overlay of Macrophages (Macs) and

Figure 6 :

 6 Figure 6: Each tumor archetype is defined by a unique combination of immune gene expression pattern A. Bubble plot of NK cells (natural killer cells), B cells, plasma cells and Mast cells associated gene expression in the live compartment grouped by clusters/archetypes identified in the IPI cohort using 10-feature clustering. B. Heatmap and hierarchical clustering of the median chemokine gene expression of all chemokines in the Chemokine phenotype signature, grouped by cluster/archetype. C. Bubble plot of gene expression in Macrophages (M1, M2) and Dendritic cells (Co-Stim, DC) function in the myeloid compartment, grouped by cluster/archetypes in the IPI cohort. D. Heatmap and hierarchical clustering of the median gene expression of B cells (all viable RNAseq compartment), NK cells (all viable RNAseq compartment), plasma cells (all viable RNAseq compartment) T cells phenotypes (Tconv

Figure 7 :

 7 Figure 7: Immune archetypes tie closely to tumor biology and disease outcome. A. A. Left-UMAP display and graph-based clustering of tumor immune archetypes using 10-feature clustering color-coded by tumor type. Each dot represents a single patient. Right-stacked bar plot of the tumor type distribution for 10-feature immune archetypes. B. UMAP overlay (Left) and box and whisker plot (Right) of tumor proliferation measured by frequency of Ki67+ CD45-cells by flow cytometry in the IPI cohort. C. Bubble plot of the median gene expression, in the tumor compartment, of gene sets associated with previously identified tumor transcriptional programs grouped by cluster/archetypes in the IPI cohort. D,E. Heatmap and hierarchical clustering of immune archetype gene signatures median expression in the tumor compartment (D) and the live compartment (E) in the IPI cohort.

  Correlation plots of Tcell, Myeloid and CD90+ CD44+ Stroma gene signature scores, for each tumor specimens, against their corresponding flow score (See Fig Supp 1 and STAR protocol) color-coded according to the tumor types shown in 1C. E. Cross-whisker plots comparing median Tcell, Myeloid and CD90+ CD44+ Stroma gene signature scores by tumor type to the median flow score, color-coded according to the tumor types shown in 1C, with the interquartile range on both axes.

table S2 for details by cancer type). D.-Left-Gene score calculation method (see STAR methods).-Right-

Cell. Author manuscript; available in PMC 2023 January 06.

ACKNOWLEDGEMENTS

We thank all members of the Krummel Lab, UCSF ImmunoX, UCSF CoLabs and the UCSF Immunoprofiler Consortium for discussion and guidance while developing this study. We would like to thank Dr Nicholas Kuhn for scientific discussion. We would like to thank Dr Kenneth Hu for editing the manuscript. We would like to particularly thank Isabelle Tingin, Garry Shumakher, Elizabeth Edmiston and Meghan Zubradt for their constant support and help during this study. We would like to thank Dr Ana Catharina Silva for her help in designing the graphical abstract. Acquisition and analysis of certain human samples described in this study was partially funded by contributions from AbbVie, Amgen, Bristol-Myers Squibb, and Pfizer as part of the UCSF Immunoprofiler Initiative. Further support came from the NIH (R01CA197363 and U01CA217864) for MFK, NCI R01 CA178015, CA222862, CA227807, CA239604, CA230263], U24 [CA210974], U54 [CA224081] for EAC and P30 : P30CA082103 for ABO. Finally, we thank all patients and their families, for placing their trust in us.

CONSORTIA

The UCSF Immunoprofiler Initiative (IPI) includes: Matthew Spitzer, Lawrence Fong, Amanda Nelson, Raj Kumar, Justin Lee, Arun Burra, Joy Hsu, Caroline Hackett, Karen Tolentino, Jasmine Sjarif, Peter Johnson, Evans Shao, Darrell Abrau, Leonard Lupin, Cole Shaw, Zachary Collins, Tasha Lea, Carlos Corvera, Eric Nakakura, Julia Carnevale, Michael Alvarado, Kimberley Loo, Lawrence Chen, Melissa Chow, Jennifer Grandis, Will Ryan, Ivan El-Sayed, David Jablons, Gavitt Woodard, Hideho Okada, Margaret Tempero, Andrew Ko, Kim Kirkwood, Scott Vandenberg, Denise Guevarra, Erica Oropeza, Chris Cyr, Pat Glenn, Jennifer Bolen, Amanda Morton, Walter Eckalbar. Affiliations for Consortia members can be found in supplementary table S7.

Author Manuscript