Resource A single-cell tumor immune atlas for precision oncology Paula Nieto,1 Marc Elosua-Bayes,1 Juan L. Trincado,1 Domenica Marchese,1 Ramon Massoni-Badosa,1 Maria Salvany,2 Ana Henriques,2 Juan Nieto,1 Sergio Aguilar-Fernández,1 Elisabetta Mereu,1 Catia Moutinho,1 Sara Ruiz,1 Patricia Lorden,1 Vanessa T. Chin,3,4,5 Dominik Kaczorowski,3 Chia-Ling Chan,3 Richard Gallagher,5,6 Angela Chou,7,8,9 Ester Planas-Rigol,10 Carlota Rubio-Perez,10 Ivo Gut,1 Josep M. Piulats,11 Joan Seoane,10,12,13,14 Joseph E. Powell,3,15 Eduard Batlle,2,12,14 and Holger Heyn1,16 1 CNAG-CRG, Centre for Genomic Regulation (CRG), Barcelona Institute of Science and Technology (BIST), 08028 Barcelona, Spain; 2 Institute for Research in Biomedicine (IRB Barcelona), Barcelona Institute of Science and Technology, 08028 Barcelona, Spain; 3 Garvan-Weizmann Centre for Cellular Genomics, Garvan Institute of Medical Research, Darlinghurst NSW 2010, Sydney, Australia; 4 St. Vincent’s Hospital Clinical School, University of New South Wales, Sydney, Australia, and St. Vincent’s Hospital Sydney, Darlinghurst NSW 2010, Australia; 5St. Vincent’s Hospital Sydney, Darlinghurst NSW 2010, Australia; 6University of Notre Dame, Chippendale NSW 2007, Sydney, Australia; 7Cancer Diagnosis and Pathology Group, Kolling Institute of Medical Research, Royal North Shore Hospital, St. Leonards NSW 2065, Australia; 8NSW Health Pathology, Department of Anatomical Pathology, Royal North Shore Hospital, Sydney NSW 2065, Australia; 9University of Sydney, Sydney NSW 2006, Australia; 10Vall d’Hebron Institute of Oncology (VHIO), Vall d’Hebron University Hospital, 08035 Barcelona, Spain; 11Medical Oncology Department, Institut Català d’Oncologia - ICO; Clinical Research in Solid Tumors Group - CREST, Bellvitge Biomedical Research Institute IDIBELL-OncoBell; CIBERONC; 08908 L’Hospitalet de Llobregat, Barcelona, Spain; 12Institució Catalana de Recerca i Estudis Avançats (ICREA), 08010 Barcelona, Spain; 13Universitat Autònoma de Barcelona (UAB), 08193 Barcelona, Spain; 14CIBERONC, 08908 Barcelona, Spain; 15 UNSW Cellular Genomics Futures Institute, University of New South Wales, Sydney NSW 2052, Australia; 16Universitat Pompeu Fabra (UPF), 08002 Barcelona, Spain The tumor immune microenvironment is a main contributor to cancer progression and a promising therapeutic target for oncology. However, immune microenvironments vary profoundly between patients, and biomarkers for prognosis and treatment response lack precision. A comprehensive compendium of tumor immune cells is required to pinpoint predictive cellular states and their spatial localization. We generated a single-cell tumor immune atlas, jointly analyzing published data sets of >500,000 cells from 217 patients and 13 cancer types, providing the basis for a patient stratification based on immune cell compositions. Projecting immune cells from external tumors onto the atlas facilitated an automated cell annotation sys- tem. To enable in situ mapping of immune populations for digital pathology, we applied SPOTlight, combining single-cell and spatial transcriptomics data and identifying colocalization patterns of immune, stromal, and cancer cells in tumor sec- tions. We expect the tumor immune cell atlas, together with our versatile toolbox for precision oncology, to advance cur- rently applied stratification approaches for prognosis and immunotherapy. [Supplemental material is available for this article.] Single-cell RNA sequencing (scRNA-seq) techniques are powerful 2019; Park et al. 2020) and aging (Salzer et al. 2018; The Tabula tools for the unbiased charting of cellular phenotypes (Lafzi Muris Consortium 2020). Diseased tissues display an additional et al. 2018). Analyzing transcriptome profiles of individual cells layer of complexity, presenting cell type composition shifts and enables the fine-grained annotation of cell types and cellular newly emerging disease-specific phenotypes (Ramachandran states, as well as charting the composition of complex samples. et al. 2019; Vieira Braga et al. 2019; Chua et al. 2020). In cancer, Interrogating healthy tissues at single-cell resolution provides a in addition to diverse neoplastic cell states (Patel et al. 2014; reference atlas of normal tissue organization (Regev et al. 2017) Tirosh et al. 2016b), the remodeling of the host tissue microenvi- and defines variability across individuals (van der Wijst et al. ronment has been characterized using scRNA-seq (Tirosh et al. 2018) or during development (Asp et al. 2019; Popescu et al. 2016a; Puram et al. 2017). Phenotyping single cells from the tumor microenvironment (TME) has led to the identification of cancer- Corresponding author: holger.heyn@cnag.crg.eu Article published online before print. Article, supplemental material, and publi- © 2021 Nieto et al. This article, published in Genome Research, is available un- cation date are at https://www.genome.org/cgi/doi/10.1101/gr.273300.120. der a Creative Commons License (Attribution-NonCommercial 4.0 Internation- Freely available online through the Genome Research Open Access option. al), as described at http://creativecommons.org/licenses/by-nc/4.0/. 31:1913–1926 Published by Cold Spring Harbor Laboratory Press; ISSN 1088-9051/21; www.genome.org Genome Research 1913 www.genome.org Nieto et al. specific stromal cell states and supported their contribution to and annotation systems, resulting in data sets characterized by tumor progression. Functional and integrative analysis further varying granularity and nomenclature. Cell annotation is especial- support dependencies of stromal and cancer cells and their predic- ly challenging when describing novel phenotypes (e.g., cancer- tive value for patient outcome. In particular, cancer-associated fi- specific cell states), leading to discrepancies in cell labels between broblasts (CAFs) (Merlos-Suárez et al. 2011; Calon et al. 2015) studies. In this work, we generate a consensus atlas of the tumor and tumor-resident immune cells (Fridman et al. 2012; Yofe et al. immune microenvironment through the integration of published 2020) have been identified as biomarkers for patient stratification scRNA-seq data sets from 13 cancer types for a harmonized and au- and, importantly, as an effective target for therapeutic interven- tomated cell annotation, as an input for patient stratification and tion in oncology (Tumeh et al. 2014; Tauriello et al. 2018). to spatially localize immune cell populations in tumor sections. Clinically most advanced, immune cells are now the target of im- munotherapy (e.g., immune checkpoint inhibitors, ICI), stimulat- ing the immune system to recognize and destroy cancer cells. Results Single-cell transcriptomics has critically advanced our under- standing of immune cell diversity in the TME by generating high- Generating a tumor immune cell atlas resolution landscapes of different cancer types. Combined with T Single-cell transcriptome profiling of tumors provides an unbiased cell receptor (TCR) genotyping and receptor-ligand screening overview of the heterogeneity of cancer cells and their microenvi- (Mimitou et al. 2019), scRNA-seq identified transient immune ronment. Following sample dissociation, single cells enter scRNA- cell states and dynamic tissue remodelling. Seminal studies in- seq processes either directly or following the enrichment of specific clude the scRNA-seq-based immunophenotyping of breast cancer, cell types (e.g., tumor or PTPRC+ immune cells). To generate a com- describing a tumor-specific heterogeneity and expansion of T cell prehensive tumor immune cell atlas of human cancers, we collect- states, supporting a continuous cell activation toward terminal ex- ed scRNA-seq data sets from 13 different cancer types, 217 patients, haustion, rather than discrete cellular states (Azizi et al. 2018). In and 526,261 cells. In detail, we processed data from breast carcino- melanoma, similar dynamics were shown to activate transitional mas (BC) (Azizi et al. 2018), basal cell and squamous cell carcino- and dysfunctional T cell states (Li et al. 2019). Moreover, single- mas (BCC) (Yost et al. 2019), endometrial adeno- (EA) and renal cell immunophenotyping identified tumor-specific T cell states, cell carcinomas (RCC) (Wu et al. 2020), intrahepatic cholangio- which are predictive for ICI-based therapy outcome (Sade- (ICC) and hepatocellular carcinomas (HCC) (Ma et al. 2019; Zhang Feldman et al. 2018). A lung cancer study identified an immune ac- et al. 2019), colorectal cancers (CRC) (Lee et al. 2020; Wu et al. tivation module characterized through high frequencies of 2020), pancreatic ductal adenocarcinomas (PDAC) (Peng et al. PDCD1+ CXCL13+ activated T cells, IgG+ plasma cells, and SPP1+ 2019), ovarian cancers (OC) (Schelker et al. 2017), non-small-cell macrophages, introducing the concept of an immune cell compo- lung cancers (NSCLC) (Lavin et al. 2017; Lambrechts et al. 2018; sition-based patient stratification to complement current genomic Wu et al. 2020), and cutaneous (CM) and uveal (UM) melanomas (e.g., mutational load) or biomarker (e.g., CD274 or CTLA4) strat- (Supplemental Table 1; Fig. 1A; Sade-Feldman et al. 2018; Li et al. egies (Leader et al. 2020). 2019; Durante et al. 2020). For a cell type stratification and consis- In more general terms, single-cell sequencing-based immuno- tent annotation, the data sets were analyzed separately before join- phenotyping identified cancer-specific states and composition bi- ing immune cells into a pan-cancer tumor immune reference atlas ases across all major immune cell types that colocalize with cancer (Supplemental Figs. 1, 2). In line with previous findings, immune cells (Qian et al. 2020). Major alterations have been described for T cells clustered by cell identity rather than patient origin, allowing cells (Yost et al. 2019), B cells (Helmink et al. 2020), tumor-associ- the straightforward subsetting of the immune compartment. ated macrophages (TAMs) (Lee et al. 2020), and dendritic cells We hypothesize that joining cells from different cancer types (Lavin et al. 2017), suggesting a global perturbation of the immune into a single reference data set may define commonalities and har- system in cancer. Comprehensively understanding the causes and monize annotations between studies. To this end, we integrated consequences of perturbed immune cell function across cancer 317,111 immune cells using canonical correlation analysis (Fig. types could provide the basis to identify novel therapeutic targets 1B; Butler et al. 2018). This approach identifies common cellular and could lay the groundwork for an immune-based patient strat- phenotypes and allows merging data sets from different studies ification in precision oncology. Current ICI stratification strategies and technologies for joint analyses. Data integration with alterna- involve the assessment of the tumor mutational burden or inflam- tive tools resulted in comparable results, though with reduced iLISI mation signatures. However, despite having favorable immune integration score (Korsunsky et al. 2019) and data set-specific tech- profiles, many tumors do not respond to treatment, suggesting ad- nical clusters (Supplemental Fig. 3). Following integration, cells ditional mechanisms that confer resistance to therapy. In this re- separated into 25 clusters representing major immune cell types, gard, the spatial distribution of immune cells has proven to be including 12 T cell, five macrophage/monocyte, three dendritic important for ICI response, with excluded tumors blocking effec- cell (DC), three B- and plasma B cell, one natural killer (NK), and tive immune cell action, despite the presence of favorable cell one Mast cell cluster (Fig. 1B,C). To test the robustness of the clus- types at their boundaries (Chen and Mellman 2017). In contrast, ters and their associated signatures, we trained a random forest (RF) tumors invaded by reactive immune cells that clonally expand classifier to predict cell annotation based on the 25 signatures and have shown increased response rates. However, current spatial im- performed a fivefold cross-validation to assess biases and variance. mune profiling approaches enable only targeted profiling of The mean accuracy and kappa statistic across folds were 0.76 (sd = mRNAs or proteins. On the other hand, the latest spatial transcrip- 0.0048) and 0.75 (sd = 0.0050), respectively, a 3.0-fold and 4.2-fold tomics (ST) techniques (Ståhl et al. 2016; Rodriques et al. 2019) increase with respect to random signatures and comparable to val- provide unbiased transcriptome-wide profiles but average gene ex- ues obtained in other high-quality atlases (Supplemental Fig. 4A,B; pression profiles of multiple cells. Zeisel et al. 2018). In addition, both the interquartile range (IQR) More broadly, it remains challenging to generalize findings and the range of accuracies were lower when using cell type–specif- across cancer types. There are no standardized analysis pipelines ic signatures, suggesting a low variance across test sets. Finally, 1914 Genome Research www.genome.org Tumor immune atlas for precision oncology A B C E F G D Figure 1. Characterization of the tumor immune cell atlas. (A) Number of cells (top) and patients (bottom) per cancer type included in the atlas. (B) UMAP of 317,111 immune cells from 13 cancer types colored by annotated cell type. (C) Total number of cells of each immune cell type/state; color code as in A. (D) Marker gene expression levels for broader cell types (left) and only T cells states (right); color code as in A. (E) Cancer type proportions for each annotated cell type/state. (F ) Number of unique patients representing each cell type/state in the atlas. (G) Expression of the top four differentially expressed genes per cell type/state; colored by cell type (as in A) and cancer type (as in E). Genome Research 1915 www.genome.org Nieto et al. stratifying the accuracies by cell type revealed that the low degree larity assessment across the 13 cancer types (Fig. 2A). Samples of of misclassification corresponded to cell types within the same lin- different tumor types intermixed, and cancer type was not a eage, which share several markers (Supplemental Fig. 4C). Taken main source of variance as shown in a t-distributed stochastic together this suggests that the clustering of the atlas displays an neighbor embedding (t-SNE) of immune cell type proportions optimal balance between granularity and robustness. Of note, (Fig. 2B). In contrast, abundance of CD4 effector memory T cells deeper clustering (31 clusters) resulted in a more fine-grained an- explained most variance in the data set (principal component 1, notation of mainly the B cells, with subpopulations representing 43.7%), followed by the relative frequencies of C1QC TAMs and naive, activated, memory, unswitched memory, and proliferative terminally exhausted CD8 T cells (32.6% and 9.5%, respectively) subtypes (Supplemental Fig. 5). (Fig. 2C). Clusters were manually annotated using canonical markers Correlating immune cell population frequencies identified an and curated gene signatures that defined their identities (Supple- intriguing pattern of tumors being mutually exclusive for the pres- mental Tables 2, 3). Exemplarily, CD4 T cells representing regula- ence of lymphoid or myeloid cell types (Supplemental Figs. 7, 8A). tory, helper, and naive states were defined by FOXP3, CXCL13, and Moreover, lymphoid-enriched tumors presented either naive- SELL marker gene expression, respectively (Fig. 1D; Supplemental memory, effector-memory and recently activated profiles, or Fig. 6). Cytotoxic CD8 T cells expressed high levels of GZMA, more differentiated (regulatory, cytotoxic, and exhausted) pheno- whereas their pre- and terminally exhausted subtypes expressed types. Myeloid-enriched tumors further split into pro-inflammato- IFIT3 and HAVCR2, respectively. Macrophages split in different ry and inhibitory (C1QC and SPP1 TAMs) subtypes, the former subtypes, such as C1QC, SPP1, and pro-inflammatory states, correlating with the presence of plasmacytoid and myeloid DCs. with C1QC, SPP1, and CXCL8 as respective marker genes. We fur- A hierarchical k-means clustering using immune cell proportions ther observed actively dividing lymphoid and myeloid cell types as features defined six clusters with largely different compositions with MKI67 and STMN1 as common markers for proliferation. Plas- (Fig. 2A,B,D,E). In spite of the different sizes of the clusters, almost ma B cells specifically expressed IGKC, which distinguished them all cancer types were presented in each cluster, confirming large from naive/memory B cells expressing MS4A1. The large number commonalities of immune cell compositions between them (Fig. of immune cells enabled the identification of rare cell states, 2F). Tumors in cluster 1 (C1) showed a high proportion of recently such as proliferative B cells (834 cells) and macrophages (1048 activated and naive-memory CD4 T cells, whereas C2 had high cells) as well as plasmacytoid and conventional DC populations amounts of cytotoxic and terminally exhausted CD8 cells (Fig. (2140 and 2238) (Fig. 1C). All cell types were present in multiple 2D,E). C3 displayed exceptionally high levels of macrophages cancer types and patients (Fig. 1E,F), which supports the robust- (C1QC and SPP1), C4 was driven by increased frequencies of plas- ness of the integration approach. Our analyses revealed that sever- ma B cells, C5 contained high proportions of CD4 effector memo- al cell states such as naive, proliferative, transitional memory and ry T cells, and C6 was very high in B cells. To confirm the terminally exhausted T cell subtypes, and all macrophage states proportional differences of cell types between clusters, we applied were abundant in all tumors, highlighting common mechanisms scCODA, confirming the enrichment of cluster-specific cell states that operate in the TME independently of the tissue of origin and the largely mutually exclusive distribution of myeloid and and mutational background. Yet, we also found cell states enriched lymphoid cell types in tumors (Supplemental Fig. 7). The cell- in specific cancer types, suggesting cancer-specific immune envi- type enrichment was further confirmed by generalized linear mod- ronments (Fig. 1E). In particular, effector memory CD8 T cells els interrogating the differences in cell abundance across experi- were frequently found in EA and NSCLC, proliferating B cells mental conditions and cancer types, also accounting for batch were very abundant in PDAC, whereas Mast cells were mainly de- effects. We observed the cluster-defining cell types to be statisti- tected in BC. In summary, we built a catalog of immune cell types cally enriched within their respective clusters (Supplemental Table and state markers present in the TME of multiple cancer types (Fig. 5). Consistent with the distinct cell type proportions, the clusters 1G; Supplemental Table 4). This resource may help annotate future presented specific gene expression signatures (Supplemental Fig. single-cell tumor data sets. 8B; Supplemental Table 6). Intriguing from a therapeutic perspec- tive, C5 showed a striking increase of exhaustion markers on CD8 cells with significantly higher levels of LAG3, PDCD1, and CTLA4 Tumor stratification by immune cell composition (ANOVA, P < 0.01) (Supplemental Fig. 8C,D). The fact that these Tumors of the same cancer type have been described to be hetero- markers are also target for immune-therapy suggests C5 tumors geneous, presenting distinct genetic and epigenetic alterations as to be more susceptible to ICI treatment but could also shed light well as gene expression signatures, allowing their stratification on the role of CD4 effector memory cells, strongly enriched in into subtypes. Cancer subtypes have a clear impact on clinical these patients (Fig. 2D,E). In line with this, the presence of specific management, being predictive for patient prognosis and therapy immune cell states, such as regulatory T cells or anti-inflammatory response (e.g., consensus molecular subtypes in CRC [CMS]) TAMs, has been linked to immune-therapy efficacy. Thus, it is (Guinney et al. 2015). To date, gene expression subtypes are de- tempting to speculate that the here-identified patient clusters dif- fined through the analysis of bulk tumor samples and are frequent- fer in their response to ICIs. ly characterized by gene signatures of stromal or immune cell types To facilitate the classification of immune profiles of future (Calon et al. 2015). Single-cell resolved tumor maps confirmed this data sets, we trained a RF classifier with the 25 immune cell popu- contribution of the TME to the subtype classification and allowed lation achieving a highly accurate classification (accuracy: 0.8, an even more fine-grained interpretation of subtype composition 95% CI: [0.593, 0.932]; P-value: 5.36 × 10−5) (Supplemental Meth- (Lee et al. 2020). ods). Of note, the cluster-specific cell types were also the most im- Detecting varying immune cell compositions within cancer portant variables when training the RF classifier using cell type types but conserved profiles across cancers, we sought to establish proportions as features to predict cluster identities (Supplemental a pan-cancer immune classification system. We used immune cell Fig. 9). To integrate our patient clusters with large clinical RNA se- type and state frequencies of the reference atlas as input for simi- quencing cohorts, we used immune cell types defining our single- 1916 Genome Research www.genome.org Tumor immune atlas for precision oncology A B C D E F Figure 2. Patient stratification based on the tumor immune cell composition. (A) Cell type composition of patients colored by cell type/state frequencies. Patients are clustered (C1–6) into groups of similar cell type composition. Cancer and cluster identities are indicated below. (B) Dimensionality reduction representation (t-SNE distribution) of cell type frequencies in cancer patients colored by cancer type (top) and cluster identity (bottom). (C) Variance con- tribution to the two first principal components (PCs) of the top variable cell types. (D) Frequencies (% total cells) of cell types representative for cluster 1–6. (E) Heat map representation of cell type frequencies within each cluster. (F) Cancer type contribution to the six immune clusters; color code as cancer types in A. cell derived clusters (Fig. 2C) to deconvolute bulk RNA sequencing (LUAD), 472 skin cutaneous melanoma (SKCM), and 521 colon ad- data sets from the Tumor Cancer Genome Atlas (TCGA) program. enocarcinoma (COAD) samples. Following deconvolution, we ap- Specifically, we inferred tumor immune cell proportions for 1222 plied the aforementioned RF classifier to assign TCGA samples to breast invasive carcinoma (BRCA), 594 lung adenocarcinoma the six single-cell derived patient immune subtypes. Here, we Genome Research 1917 www.genome.org Nieto et al. could assign patients from all cancer types to the single-cell de- D). Again, distinct T cell states emerged from the dimensional re- rived clusters, with the exception of our smallest cluster enriched duction plots or after clustering analysis (Fig. 4A–C). Consistently, in plasma B cells (Supplemental Fig. 10). cell projection of the experimentally enriched T cell fraction to the reference atlas identified defined T cell states (Fig. 4E). Clonal ex- pansion according to TCR sequencing could only be detected in A resource for immune cell annotation cells assigned to be proliferative, cytotoxic, or exhausted, whereas naive, naive-memory, and recently activated cells were not of clon- To demonstrate the predictive value of the atlas, we generated que- al origin (Fig. 4E,F). ry data sets from different cancer types and varying experimental designs. Following clustering of the query data sets, we projected either single cells or clusters onto the atlas using a reference-based Spatial localization of immune cells in tumor sections projection (Butler et al. 2018) or cluster matching (Mereu et al. Widespread evidence indicates that the spatial distribution of im- 2020) tools, respectively. As proof-of-concept, we performed mune cells is important for ICI response (Chen and Mellman scRNA-seq for human primary OC and UM, liver metastases 2017). To explore this issue, we combine single-cell reference atlas (from primary UM), as well as a brain metastasis (from primary immune profiles with ST data from tumor sections. The aim of this CM). approach is to provide spatial tumor maps that precisely delineate Using cell-by-cell projection, we matched UM and OC cells to colocalization of immune cells in tumors. To integrate both data specific cell types and states of the atlas reference (Fig. 3A,D). modalities, we applied SPOTlight, a nonnegative matrix factoriza- Query cells projected to macrophage and T cells, with a large vari- tion (NMF)-based spatial deconvolution framework (Elosua-Bayes ety of cell states being detected. The UM harbored all macrophage et al. 2021). SPOTlight identifies cell type–specific topic profiles subtypes, including abundant SPP1 and C1QC TAMs and a small (gene expression signatures) from scRNA-seq data, which are sub- fraction of proliferative macrophages (Fig. 3A). In addition, the en- sequently used to deconvolute ST spots. To predict and quantify tire spectrum of T cell phenotypes could be assigned, with a high the location of the 25 cell types and states of the tumor immune fraction of cells projecting to the terminally exhausted CD8 T reference, gene expression profiles were translated into SPOTlight cell state. The OC data set contained less SPP1 TAMs and exhibited topic profiles. Topics were highly specific to clusters, resulting in an increased proportion of proliferative macrophages (Fig. 3D). We high specificities to localize immune cell types (Fig. 5A). To avoid also detected a small number of conventional DC and monocytes. conflicting signals and inaccurate predictions related to the pres- We next used clustering analyses on the query data sets and subse- ence of strong cell cycle signatures from cancer cells, we removed quently matched the identified clusters to the atlas reference. proliferating immune cells from the deconvolution step and There were seven and three clusters for UM and OC, respectively included a general proliferation signature to identify areas of divid- (Fig. 3B,E). Matching the query and reference atlas clusters resulted ing cells. To map immune cells within tumor sections, we generat- in a clear separation between the lymphoid and myeloid popula- ed ST data sets for an oropharyngeal squamous cell carcinoma tions (Fig. 3C,F). In UM, cytotoxic and terminally exhausted cells (SCC) metastasis and analyzed two sequential sections of an inva- were enriched in cluster 1 and cluster 6, in line with their distribu- sive ductal BC (Visium platform, 10x Genomics). tion in the dimensionality reduction representation (Fig. 3A). SPP1 Clustering the spatial SCC data identified four distinct tran- TAMs showed higher scores in cluster 2, whereas C1QC TAMs fell scriptional and proliferative active areas containing cancer cells into cluster 4, an assignment consistent with the cell-by-cell pro- (cluster 1/2), surrounded by stroma (cluster 0) and an area enriched jection. Cluster 7 matched to both lineages, likely depicting cell in immune cells (cluster 3) (Fig. 5B–E). Mapping the immune cell doublets. distribution in the sections using SPOTlight returned proportions The projection of cells from UM liver metastases rendered a for each cell within a spot, predictions that were validated using an precise separation of T cell states (Fig. 3G), which split into naive, alternative spot deconvolution tool (Supplemental Fig. 11; Cable recently activated, cytotoxic, and proliferative subpopulations. We et al. 2021). This resulted in a clear regionalization of distinct im- found a distinct frequency of exhausted cells between both UM mune cell types and states (Fig. 5D; Supplemental Fig. 12). The metastases. Although much less abundant, all DC and B cell states cancer areas cluster 1/2 presented similar patterns, with an enrich- could also be assigned in both samples. Cell-by-cell projection of T ment of the proliferation signature and SPP1 macrophages (Fig. 5F, cells from a NSCLC brain metastasis also identified most cell states, G; Supplemental Fig. 13). Cluster 3, in turn, presented a specific with a high proportion of terminally exhausted CD8 T cells (Fig. immune infiltration pattern dominated by an enriched presence 3H). Combined scRNA-seq and TCR genotyping assigned clonally of B cells and cDCs. The stroma section, cluster 0, harbored regula- expanded T cells as being either actively proliferating or terminally tory T cells and terminally exhausted CD8 T cells and was specifi- exhausted, in line with the expected history of tumor reactivity cally enriched in C1QC TAMs. An interaction matrix computed and proliferation of exhausted cells (Fig. 3H,I; Azizi et al. 2018; based on the predicted proportions on the spatial distribution of Li et al. 2019). cell states in the section highlighted colocalization and mutual ex- We further wondered about the applicability of the atlas as a clusive patterns (Fig. 5H; Supplemental Fig. 14). Here, spatial reference across species. To tackle this question, we generated colocalization provided further insights into immune-modulating scRNA-seq data sets for two liver metastases derived from mouse mechanisms of cancer and immune-suppressive stromal cells (e.g., CRC organoids (Tauriello et al. 2018), one PTPRC selected and regulatory T cells). Exemplarily, recently activated T cells colocal- the other enriched in T cells, and projected individual cells onto ized with effector/naive memory T cells, whereas they were exclud- the human reference. The mouse data sets could be projected ed from areas with regulatory and exhausted T cells as well as with high confidence (∼70% of cells with matching probabilities C1QC TAMs and plasma B cells. Accordingly, regulatory T cells > 0.5). Mouse clusters of main subtypes could be readily annotated, showed a positive correlation with cytotoxic cells (including and specific subpopulations could also be assigned to distinct T/B NK), pre- and terminally exhausted CD8 T cells, and C1QC cell or macrophage cell states using the human reference (Fig. 4A– TAMs, but a negative correlation with recently activated and 1918 Genome Research www.genome.org Tumor immune atlas for precision oncology A B C D E F G H I Figure 3. Automated annotation of external human tumor-derived immune cells using the tumor immune atlas as reference. (A,B) UMAP representation of immune cell transcriptomes from a primary uveal melanoma colored by their predicted cell type/state based on the tumor immune reference (A) or using unsupervised clustering (B). (C) Marker correspondence (Jaccard index) between uveal melanoma clusters (B) and the cell type clusters of the reference atlas. (D,E) UMAP representation of immune cell transcriptomes from a primary ovarian carcinoma colored by predicted cell type/state (D, color code as in A) and after clustering (E). (F ) Marker correspondence (Jaccard index) between the ovarian cancer clusters and the cell type clusters of the reference atlas. (G) UMAP representation of immune cells from two uveal melanoma liver metastasis colored by their predicted cell type (color code as in A). (H,I ) UMAP representation of T cells isolated from a brain metastasis colored by their predicted cell type (H, color code as in A) and clonal expansion profiled through TCR genotyping (I ). Genome Research 1919 www.genome.org Nieto et al. A B C D E F Figure 4. Projecting mouse tumor immune cells onto the human reference atlas. (A,B) UMAP single-cell transcriptome representations of mouse T cells isolated from tumor organoids. Cells are color-coded by their predicted cell type based on the human reference atlas (A) or by cluster identity (B). (C ) Cell type composition of each cluster (color code as in A). (D) Marker gene correspondence (Jaccard index) between the mouse immune clusters and the cell type clusters of the reference atlas. (E,F ) UMAP representation of mouse T cells isolated from tumor organoids colored by predicted cell type (E, color code as in A) and clonality based on expanded TCR clonotypes (F). effector/naive memory T cells. It has been shown that regulatory T In order to validate the colocalization of abovementioned cell cells inhibit effector T cells by releasing ENTPD1 (previously type and states, we performed a digital in situ hybridization (digital known as CD39) and NT5E (also known as CD73) that bind to ISH) experiment by correlating the presence of canonical expres- the A2A receptor (Allard et al. 2020). In line with this finding, sion markers of these cells. Although being less sensitive, ST gener- and in addition to the exclusive pattern of regulatory and activated ates accurate molecule counts, thus returning an adequate signal and effector T cells, we detected a positive expression correlation of for targeted colocalization (same spot) analysis. To enhance the ENTPD1 and NT5E with the regulatory T cell predicted proportions signal, we denoised the data (van Dijk et al. 2018) before comput- (P < 2.2 × 10−16 and P = 2.5 × 10−6, respectively). Along with the de- ing a clustered correlation matrix of marker genes (Supplemental tected spatial organization with other immune cell types, regulato- Fig. 18A). We observed strong correspondence of cell type marker ry T cells illustrate the value of spatial immune mapping to gene colocalization with the predicted proximity of respective cell advance our understanding of tumor immunology. Analysis of a types. Specifically, the localization of GNLY and NKG7 (cytotoxic T second SCC case confirmed the regional restriction of tumor im- cells) correlated significantly with FOXP3 and CTLA4 (regulatory T mune cells and confirmed the previously described colocalization cells) (Supplemental Fig. 18B–D). On the other hand, the correla- and mutually exclusive patterns (Supplemental Figs. 15–17). tion between SPP1, CLEC7A (SPP1 TAMs), and cancer (EPCAM), 1920 Genome Research www.genome.org Tumor immune atlas for precision oncology A B D C F E G H Figure 5. Spatial mapping of the reference immune cell types using ST oropharyngeal squamous cell carcinoma (SCC) sections. (A) Cell type–specific topic profiles presenting a high topic/cell type specificity. (B) ST profiled section of a SCC primary tumor. Tissue stratification according to unsupervised clustering. (C ) The number of unique molecular identifiers (UMI) recovered from each spot indicate the areas transcriptionally most active. (D) Pie chart representation showing proportions (per ST spot) of SPOTlight-predicted immune cells based on the single-cell immune reference atlas. To visualize spa- tially variable cell types, only immune cell types present in <75% of the spots are displayed. (E) UMAP embedding of ST spots presenting the cell cycle phase (left) and cluster identity (right) for each spot. (F ) Box plots displaying significant differences in cell type proportion of clusters (ANOVA test). (G) Location and proportion of significantly differentially located cell types in the SCC section. (H) Clustered correlation matrix between the predicted cell type propor- tions identifying colocalization (red) and exclusive (blue) immune distribution patterns. Genome Research 1921 www.genome.org Nieto et al. as well as proliferation markers (STMN1, MKI67, CDK1) confirmed cells (PDCD1, LAG3, HAVCR2), and C1QC macrophages (MRC1, the colocalization of SPP1 TAMs with proliferating cancer cells CD163) (Supplemental Fig. 18A). (Supplemental Fig. 18E–G). Finally, our digital ISH experiment A cancer-specific regional distribution was also observed in confirmed the colocalization of stromal cells (ACTA2, PDGFRA) the invasive ductal BC and replicated over serial sections of the tu- with regulatory T cells (FOXP3, CTLA4), terminally exhausted T mor (Fig. 6A–C; Supplemental Fig. 19). The gene expression A E B F C D Figure 6. Tumor immune reference mapping ST section from a ductal breast carcinoma. (A) Pie chart representation showing proportions (per ST spot) of SPOTlight-predicted immune cells based on the single-cell immune reference atlas. To visualize spatially variable cell types, only immune cell types pre- sent in <75% of the spots are displayed. (B,C) Estrogen receptor 1 (ESR1, B) and erb-b2 receptor tyrosine kinase 2 (ERBB2, also known as HER2, C) gene expression levels on the ST section, indicating profound regionality of the expression. (D) Tissue stratification and labeling according to unsupervised clus- tering. (E) Box plots of significantly differentially localized cell type proportions between the clusters (ANOVA test). Differences between tumor areas (i.e., HER2 + and ESR1 +) are observed, suggesting differential tumor microenvironments of the tumor subclones. (F) Location and proportion of immune cell types with local enrichment. 1922 Genome Research www.genome.org Tumor immune atlas for precision oncology profiles of cancer cells pointed to a strong subclonal structure, with mor immune cell reference, in combination with SPOTlight- mutually exclusive localizations of BC subtype-specific ERBB2 based spatial mapping of ST data, enables harmonization of the (HER2+) and ESR1 (ESR1 +) positive clones. This subclonal architec- spatial assessment of tumor infiltrating immune cells, a first ture was directly associated with local enrichment of distinct cell step toward an automated digital pathology framework. We fore- states (Fig. 6D; Supplemental Figs. 20, 21). The ESR1 +/HER2 − see the regional distribution of immune cell types becoming an region showed the increased presence of proliferating cells and important feature for the prediction of immunotherapy outcome. SPP1 macrophages (Fig. 6E,F; Supplemental Figs. 22, 23). In con- Moreover, longitudinal sampling for spatial immune cell map- trast, the HER2 +/ESR1 − region showed a higher proportion of reg- ping throughout ICI treatment could aid in explaining therapy ulatory and cytotoxic T cells, C1QC TAMs, and NK cells. In this action and point to further weaknesses approachable through region, coupled to the higher proportion of regulatory T cells, we combinatorial therapy strategies. However, it is of note that pre- also observed an enrichment of pre-exhausted and terminally ex- dicting the presence of immune cell types from a reference atlas hausted CD8 T cells, different inactivation stages of the immune cannot comprehensively take into account cell type– or patient- cell response. Lastly, we observed an enrichment of plasma B cells specific peculiarities. Overlapping gene expression signatures, as and C1QC TAMs lining the cancer areas in the fibrotic tissue. The shown for the proliferation-related transcriptional program, can importance of a quantitative assessment of cell types within ST challenge distinguishing the presence of immune cells from cell spots is displayed by the ubiquitous presence of SPP1 macrophages types that are not present in the training set. Consequently, ad- across all areas but the significant proportional enrichment in spe- ditional analysis, such as the here-applied digital ISH validation, cific tumor fractions (Fig. 6E,F). Similar to the SCC, spatial correla- or experimental validation should be considered. tion analysis in BC, among other features, identified the In conclusion, our analysis allowed us to catalog tumor-resi- colocalization of proliferating cells and SPP1 macrophages along dent immune cells across a large patient cohort and to draw a com- with the anticorrelation with regulatory and cytotoxic cells and prehensive tumor immune cell atlas of human cancers. We C1QC TAMs (Supplemental Fig. 14B,C). envision the application of the atlas in multiple areas of precision oncology. Therefore, we provide a framework for an immune- based patient stratification, the feasibility to integrate newly gen- Discussion erated patient single-cell data, and a toolbox to map immune cells We jointly analyzed 317,111 immune cells from 217 patients, in- directly in tumor sections. Following the integration of clinical co- cluding major (e.g., T cells, B cells, and macrophages) and rare (e.g., hort single-cell studies with patient outcome and response meta- proliferating and dendritic cells) cell types and states infiltrating data, we expect the atlas to be predictive for patient prognosis human tumors. We stratified cell types into subpopulations and and immunotherapy response to a level that greatly exceeds cur- annotated a fine-grained map of immune cell states. We envision rently applied stratification strategies. this atlas to serve as a reference for a harmonized annotation of future data sets, which we tested by projecting single cells from pri- mary human tumors and metastases as well as mouse immuno- Methods therapy models onto the atlas. The here-identified 25 cell types/ states present a compendium of immune cell populations shared Data preprocessing across cancer types to enable personalized medicine applications. We initially preprocessed each scRNA-seq data set independently: Whereas the RF classifier pointed to the robustness of the clusters, breast carcinomas (obtained from the NCBI Gene Expression Om- increasing the cluster number could identify additional biological nibus [GEO; https://www.ncbi.nlm.nih.gov/geo/] accession num- relevant and functional cell states. This was exemplified by the bers GSE114727 and GSE114725), basal cell and squamous cell identification of additional B cell states through deeper clustering, carcinomas (GEO GSE123814), endometrial adeno- and renal and we expect the large cell number of this atlas to allow an even cell carcinomas (GEO GSE99254), intrahepatic cholangio- and deeper annotation. hepatocellular carcinomas (GEO GSE125449, GSE140228), colo- Observing high cell state heterogeneity across patients and rectal cancers (GEO GSE132465, GSE99254), pancreatic ductal cancer types enabled a patient stratification based on the immune adenocarcinomas (Genome Sequence Archive [GSA; https://ngdc composition. This immune-based classification identified pan- .cncb.ac.cn/gsa/] accession number CRA001160), ovarian can- cancer subtypes exclusively enriched in recently activated CD4 T cers (https://figshare.com/s/711d3fb2bd3288c8483a), non-small- cells, cytotoxic/terminally exhausted CD8 T cells, CD4 effector cell lung cancers (ArrayExpress [https://www.ebi.ac.uk/array memory, C1QC tumor-associated macrophages, or B cells and express] accession numbers E-MTAB-6149, E-MTAB-6653 and plasma B cells. We believe pan-cancer immune classifications to GEO GSE97168, GSE99254), and cutaneous (GEO GSE120575, GSE123139) and uveal (GEO GSE139829) melanomas. be highly informative and clinically predictive when integrating From the original count matrices, we filtered the low-quality tumor samples from single-cell studies with patient survival and cells based on mitochondrial RNA percentage, number of UMIs response to therapy metadata. We thus envision our immune- (unique molecular identifiers), and number of different genes based patient classifier to be a useful tool for future patient stratifi- (thresholds adjusted separately for each data set). Then, we ob- cation, predicting prognosis and immunotherapy response. tained the highly variable genes, corrected, logarithmically nor- Moreover, using the classifier, the pan-cancer immune classifica- malized, and scaled the counts (default parameters) of each data tion system could be extended to additional cancer types and drive set using Seurat (version 3.2.0) (Stuart et al. 2019). Principal com- the design of basket clinical trials in which a common immune ponent analysis (PCA) was performed on all individual data sets. stratification and recruitment framework is applied across cancer The optimal number of principal components ranged between types. 25 and 40, depending on data set size. Dimensionality reduction Detecting regional immune cell enrichments in tumor was performed by applying the Uniform Manifold Approximation clones and adjacent areas, we foresee an application of the atlas and Projection (UMAP) algorithm, which also served as a two-di- in digital pathology to guide immunotherapy decisions. The tu- mensional embedding for data visualization. Cluster analysis was Genome Research 1923 www.genome.org Nieto et al. performed using the Louvain clustering algorithm, a method sep- for a better model. To evaluate the variable importance obtained arating different communities inside the data by optimizing the from the RF classifier, we assessed how the accuracy of the model modularity measure. The clusters were annotated by expression decreased when a variable is excluded from the model (mean of canonical markers of major cell types to consistently identify decrease accuracy) and how each variable affects the homogeneity and extract immune cells. Before data set integration, the major of the trees in the forest (mean decrease Gini). We then used cell types were identified, that is, T, B, myeloid, and dendritic cells, scCODA to statistically confirm proportional differences in cell and separated from the tumor cells. Note that, for some data sets, type abundance in the six patient clusters. Therefore, we included prior information for cell-level annotation was also available. all cell types/states of the atlas, and let the algorithm automatically select a reference cell type for scCODA’s model (transitional mem- ory CD4 T cells). Of note, using proliferative T cells as a reference Tumor immune cell atlas construction cell type produced similar results. We then tested each of the clus- To construct the atlas, we used the raw counts from the previously ters against the rest, setting the patient and the cluster as covari- annotated immune cells of each data set and followed Seurat’s ates. We used scCODA’s default model and algorithm and standard data integration pipeline. This method is based on the gathered the results of the estimated parameter and log2-fold identification of anchor cells between pairs of data sets, which change for each of the six contrasts (one per cluster). These results are used to harmonize the integrated data sets. To do this, canon- further confirmed and provided statistical support for the differ- ical correlation analysis (CCA) is applied, for which we set the ences in patient cluster cell type composition (Supplemental Fig. number of anchors to 3000. Prior to integration, we obtained the 7). To further assess the robustness of the patient clusters and data set-specific genes in order to remove them from the integra- to identify significantly enriched cell types, we used generali- tion anchors, thereby reducing noise. To do this, we joined all zed linear models (GLMs). Specifically, we used the R package the data sets and, with the FindAllMarkers function, calculated “emmeans” to perform the contrasts between the clusters and in- the genes specific for each data set. Then, we removed these genes cluded the batch (data set of origin), cancer type, and sequencing from the integration features. technologies, as covariates of the models (Supplemental Table 6). Integrating multiple data sets is a computationally intensive Both scCODA and GLM analysis confirm our observations. procedure. Due to the large number of atlas cells, we used the alter- native reciprocal PCA (RPCA) algorithm for anchor identification, as proposed by Stuart et al. (2019). This method is notably faster Data set projection onto the tumor immune cell atlas and performs considerably similarly to the standard CCA method. For the projection and annotation of external data sets onto the at- With this procedure, we built a tumor immune cell atlas of over las, we used Seurat’s anchor-transferring method. This algorithm 300,000 immune cells across 13 cancer types and 217 patients. uses the PCA structure of the reference atlas and projects it on After integration, we computed the clusters again to obtain a the query. Then, it finds pairs of anchor genes between data sets fine-grained resolution of cell types and states. We increased the that enable the projection. These anchors allow the transfer of an- resolution parameter in Seurat’s FindClusters up to 1.2, to retain notations (labels) between data sets. Here, we used an unsuper- more subtle differences between subpopulations. vised anchoring based on the first 30 principal components of We performed differential expression analysis (DEA) for all the data sets and the RNA assays with a LogNormalize normaliza- clusters to determine their marker genes. For the DEA, we used tion method. We obtained the cell-type assignment of each cell in the normalized RNA counts, instead of using the integrated data, the query data sets plus the corresponding prediction probability. to increase the accuracy identifying significant cluster makers. We also used matchSCore2’s function to correlate the markers of To further refine the cluster signatures, we computed the clusters query clusters with the markers of the 25 cell types/states of the at- using both the Wilcoxon rank-sum test and MAST (Finak et al. las. Therefore, we used the top 100 markers (average log-fold chan- 2015) methods and kept only genes identified with both ap- ge) for each of the clusters. proaches. We annotated the refined clusters using a curated gene signature of immune cell subtypes and compared it with the re- sults of the DEA analysis (using the matchSCore package [Mereu Spatial transcriptomics analysis et al. 2020]). Here, not only the major immune cell types but Clusters were annotated according to pathologists’ assessment, also different immune cell states were identified. their transcriptional activity, and differentially expressed genes. A subset of the single-cell atlas was used to train the SPOTlight model; we selected up to 100 cells coming from melanoma cancer Clustering and classification of patients based on immune for each cell type/state. Selecting cells from one of the cancer cell composition types allowed us to reduce data set-specific noise, which could con- In order to find similarities in patients’ immune microenviron- found the model. The gene set used to train the model was the ment across cancer types, we applied patient clustering using the union between the marker genes of the cell types along with the proportion of 25 immune cell states of their TME. For this, we fil- top 3000 variable genes. Marker genes for each cell type were oba- tered out patients with less than 500 cells and applied hierarchical tined with Seurat’s function FindAllMarkers, considering only pos- k-means (HKM) clustering (Chen et al. 2005; Alashwal et al. 2019) itive markers and setting the logFC and min.pct to 0 to include all using the cell types as variables and their proportion in the TME as genes. All markers were used to initialize the model basis, and unit values. We set the number of clusters to k = 6 (Supplemental Fig. variance normalization was carried out. Nonsmooth, nonnegative 9B) and used the agglomerative Ward 2 method with Euclidean matrix factorization was the method used to carry out the factori- distance metric. This algorithm iteratively joins instances in clus- zation. Cell types contributing <3% to the spot’s predicted compo- ters, reducing intra-cluster variance. Further, we built a RF classifier sition were considered fitting noise and were set as 0. We then used to predict the assignment of each patient to one of the six immune SPOTlight to map the atlas cell types to the spatial spots. For the subtype clusters. Therefore, we split the data set into training and digital ISH, we reduced data sparsity and enhanced gene expres- test sets (75% and 25% of the patients, respectively) and grew 100 sion signals with MAGIC (van Dijk et al. 2018), applied on the nor- trees to train the classifier, while automatically assessing predic- malized gene expression data. We ran MAGIC with default tor’s importance and calculating the proximity between patients parameters, except for knn set to 2, to ensure that the information 1924 Genome Research www.genome.org Tumor immune atlas for precision oncology used to recover the data structure is based on the most similar ST the co-authors. All authors read and approved the final spots. manuscript. Data access The Tumor Immune Cell Atlas count matrix and metadata References are deposited at Zenodo (https://zenodo.org/record/5186413# Alashwal H, El Halaby M, Crouse JJ, Abdalla A, Moustafa AA. 2019. The ap- .YRqbJC1h2v6). The analysis code is available at GitHub (https:// plication of unsupervised clustering methods to Alzheimer’s disease. github.com/Single-Cell-Genomics-Group-CNAG-CRG/Tumor-Immune- Front Comput Neurosci 13: 31. doi:10.3389/fncom.2019.00031 Allard D, Allard B, Stagg J. 2020. On the mechanism of anti-CD39 immune Cell-Atlas) and as Supplemental Code. Raw data are available as checkpoint therapy. J Immunother Cancer 8: e000186. doi:10.1136/jitc- stated in the Supplemental Methods. To explore the atlas, project 2019-000186 external data, and apply the immune classifier, we are currently Asp M, Giacomello S, Larsson L, Wu C, Fürth D, Qian X, Wärdell E, Custodio developing a user-friendly ShinyApp accessible at https://paula- J, Reimegård J, Salmén F, et al. 2019. A spatiotemporal organ-wide gene expression and cell atlas of the developing human heart. Cell 179: nieto-garcia.shinyapps.io/TICA/. Raw data for the projected sin- 1647–1660.e19. doi:10.1016/j.cell.2019.11.025 gle-cell and ST data sets generated in this study have been submit- Azizi E, Carr AJ, Plitas G, Cornish AE, Konopacki C, Prabhakaran S, Nainys J, ted to the NCBI Gene Expression Omnibus (GEO; https://www Wu K, Kiseliovas V, Setty M, et al. 2018. Single-cell map of diverse im- .ncbi.nlm.nih.gov/geo/) under accession number GSE158803. mune phenotypes in the breast tumor microenvironment. Cell 174: 1293–1308.e36. doi:10.1016/j.cell.2018.05.060 Butler A, Hoffman P, Smibert P, Papalexi E, Satija R. 2018. Integrating single- cell transcriptomic data across different conditions, technologies, and Competing interest statement species. Nat Biotechnol 36: 411–420. doi:10.1038/nbt.4096 Cable DM, Murray E, Zou LS, Goeva A, Macosko EZ, Chen F, Irizarry RA. H.H. is cofounder of Omniscope. J.S. is cofounder of Mosaic Bio- 2021. Robust decomposition of cell type mixtures in spatial transcrip- medicals, board member of Northern-Mosaic partnership, and tomics. Nat Biotechnol doi:10.1038/s41587-021-00830-w has ownership interests from Mosaic Biomedicals and Northern Bi- Calon A, Lonardo E, Berenguer-Llergo A, Espinet E, Hernando-Momblona X, Iglesias M, Sevillano M, Palomo-Ponce S, Tauriello DVF, Byrom D, ologics. J.S. received grant/research support from Mosaic Biomedi- et al. 2015. Stromal gene expression defines poor-prognosis subtypes cals, Northern Biologics, Roche/Glycart, and Hoffmann la Roche. in colorectal cancer. Nat Genet 47: 320–329. doi:10.1038/ng.3225 J.S. reports a scientific consultancy role for Merck Serono, GSK, Chen DS, Mellman I. 2017. Elements of cancer immunity and the cancer– immune set point. Nature 541: 321–330. doi:10.1038/nature21349 and Eli Lilly. All other authors declare no competing interests. Chen B, Harrison R, Pan Y, Tai PC. 2005. Novel hybrid hierarchical-K-means clustering method (H-K-means) for microarray analysis. In Proceedings of the 2005 IEEE Computational Systems Bioinformatics Conference - Acknowledgments Workshops(CSBW’05), pp. 105–108. Stanford, CA. IEEE Computer Society, Washington, DC. This publication is part of a project (BCLLATLAS) that has received Chua RL, Lukassen S, Trump S, Hennig BP, Wendisch D, Pott F, Debnath O, funding from the European Research Council (ERC) under the Eu- Thürmann L, Kurth F, Völker MT, et al. 2020. COVID-19 severity corre- ropean Union’s Horizon 2020 research and innovation pro- lates with airway epithelium–immune cell interactions identified by sin- gle-cell analysis. Nat Biotechnol 38: 970–979. doi:10.1038/s41587-020- gramme (grant agreement No. 810287). This work has received 0602-4 funding from the Ministerio de Ciencia, Innovación y Universi- Durante MA, Rodriguez DA, Kurtenbach S, Kuznetsov JN, Sanchez MI, dades (SAF2017-89109-P; AEI/FEDER, UE) and the Fundació La Decatur CL, Snyder H, Feun LG, Livingstone AS, Harbour JW. 2020. Marató de TV3 (201903-30-31-32). We further acknowledge fund- Single-cell analysis reveals new evolutionary complexity in uveal mela- noma. Nat Commun 11: 496. doi:10.1038/s41467-019-14256-1 ing from the St. Vincent’s Clinic Foundation (V.T.C.) and the Na- Elosua-Bayes M, Nieto P, Mereu E, Gut I, Heyn H. 2021. SPOTlight: seeded tional Health and Medical Research Council Investigator Grant NMF regression to deconvolute spatial transcriptomics spots with sin- (APP1175781, J.E.P.), the Fundación Asociación Española contra gle-cell transcriptomes. Nucleic Acids Res 49: e50. doi:10.1093/nar/ el Cáncer (AECC), FERO (EDM), Ramón Areces Foundation, Cellex gkab043 Finak G, McDavid A, Yajima M, Deng J, Gersuk V, Shalek AK, Slichter CK, Foundation, BBVA (CAIMI), the ISCIII, FIS (PI16/01278), Juan de la Miller HW, McElrath MJ, Prlic M, et al. 2015. MAST: a flexible statistical Cierva formación fellowship (C.R.-P.) and Sara Borrell fellowship framework for assessing transcriptional changes and characterizing het- (E.P.-R.). Core funding is from the ISCIII and the Generalitat de Ca- erogeneity in single-cell RNA sequencing data. Genome Biol 16: 278. talunya. We acknowledge support of the Spanish Ministry of Econ- doi:10.1186/s13059-015-0844-5 Fridman WH, Pagès F, Sautès-Fridman C, Galon J. 2012. The immune con- omy, Industry and Competitiveness (MEIC) to the EMBL texture in human tumours: impact on clinical outcome. Nat Rev Cancer partnership, the Centro de Excelencia Severo Ochoa, the CERCA 12: 298–306. doi:10.1038/nrc3245 Programme / Generalitat de Catalunya, the Spanish Ministry of Guinney J, Dienstmann R, Wang X, de Reyniès A, Schlicker A, Soneson C, Economy, Industry and Competitiveness (MEIC) through the Marisa L, Roepman P, Nyamundanda G, Angelino P, et al. 2015. The consensus molecular subtypes of colorectal cancer. Nat Med 21: 1350– Instituto de Salud Carlos III and the Generalitat de Catalunya 1356. doi:10.1038/nm.3967 through Departament de Salut and Departament d’Empresa i Con- Helmink BA, Reddy SM, Gao J, Zhang S, Basar R, Thakur R, Yizhak K, Sade- eixement. We also acknowledge the cofinancing by the Spanish Feldman M, Blando J, Han G, et al. 2020. B cells and tertiary lymphoid Ministry of Economy, Industry and Competitiveness (MEIC) structures promote immunotherapy response. Nature 577: 549–555. doi:10.1038/s41586-019-1922-8 with funds from the European Regional Development Fund Korsunsky I, Millard N, Fan J, Slowikowski K, Zhang F, Wei K, Baglaenko Y, (ERDF) corresponding to the 2014–2020 Smart Growth Operating Brenner M, Loh P, Raychaudhuri S. 2019. Fast, sensitive and accurate in- Program. tegration of single-cell data with Harmony. Nat Methods 16: 1289–1296. Author contributions: H.H. designed the study. P.N. and doi:10.1038/s41592-019-0619-0 Lafzi A, Moutinho C, Picelli S, Heyn H. 2018. Tutorial: guidelines for the ex- M.E.-B. performed all data analyses. J.L.T. analyzed TCR data. perimental design of single-cell RNA sequencing studies. Nat Protoc 13: E.M., R.M.-B., S.A.-F., and C.R.-P. supported the data analysis. 2742–2757. doi:10.1038/s41596-018-0073-y D.M., C.M., S.R., P.L., and E.P.-R. generated single-cell data sets. Lambrechts D, Wauters E, Boeckx B, Aibar S, Nittner D, Burton O, Bassez A, D.K. and C.C. generated ST data sets. M.S., A.H., and J.N. supported Decaluwé H, Pircher A, Van den Eynde K, et al. 2018. Phenotype mold- ing of stromal cells in the lung tumor microenvironment. Nat Med 24: data interpretation. V.T.C., R.G., I.G., J.M.P., J.S., J.E.P., and E.B. 1277–1289. doi:10.1038/s41591-018-0096-5 provided samples and technical and sequencing support. H.H., Lavin Y, Kobayashi S, Leader A, Amir ED, Elefant N, Bigenwald C, Remark R, P.N., and M.E.-B. wrote the manuscript with contributions from Sweeney R, Becker CD, Levine JH, et al. 2017. Innate immune landscape Genome Research 1925 www.genome.org Nieto et al. in early lung adenocarcinoma by paired single-cell analyses. Cell 169: Salzer MC, Lafzi A, Berenguer-Llergo A, Youssif C, Castellanos A, Solanas G, 750–765.e17. doi:10.1016/j.cell.2017.04.014 Peixoto FO, Attolini CS-O, Prats N, Aguilera M, et al. 2018. Identity noise Leader AM, Grout JA, Chang C, Maier B, Tabachnikova A, Walker L, Lansky and adipogenic traits characterize dermal fibroblast aging. Cell 175: A, LeBerichel J, Malissen N, Davila M, et al. 2020. CITEseq analysis of 1575–1590.e22. doi:10.1016/j.cell.2018.10.012 non-small-cell lung cancer lesions reveals an axis of immune cell activa- Schelker M, Feau S, Du J, Ranu N, Klipp E, MacBeath G, Schoeberl B, Raue A. tion associated with tumor antigen load and TP53 mutations. bioRxiv 2017. Estimation of immune cell content in tumour tissue using single- doi:10.1101/2020.07.16.207605 cell RNA-seq data. Nat Commun 8: 2032. doi:10.1038/s41467-017- Lee H-O, Hong Y, Etlioglu HE, Cho YB, Pomella V, Van den Bosch B, 02289-3 Vanhecke J, Verbandt S, Hong H, Min J-W, et al. 2020. Lineage-depen- Ståhl PL, Salmén F, Vickovic S, Lundmark A, Navarro JF, Magnusson J, dent gene expression programs influence the immune landscape of co- Giacomello S, Asp M, Westholm JO, Huss M, et al. 2016. Visualization lorectal cancer. Nat Genet 52: 594–603. doi:10.1038/s41588-020-0636-z and analysis of gene expression in tissue sections by spatial transcrip- Li H, van der Leun AM, Yofe I, Lubling Y, Gelbard-Solodkin D, van Akkooi tomics. Science 353: 78–82. doi:10.1126/science.aaf2403 ACJ, van den Braber M, Rozeman EA, Haanen JBAG, Blank CU, et al. Stuart T, Butler A, Hoffman P, Hafemeister C, Papalexi E, Mauck WM, Hao Y, 2019. Dysfunctional CD8 T cells form a proliferative, dynamically regu- Stoeckius M, Smibert P, Satija R. 2019. Comprehensive integration of lated compartment within human melanoma. Cell 176: 775–789.e18. single-cell data. Cell 177: 1888–1902.e21. doi:10.1016/j.cell.2019.05 doi:10.1016/j.cell.2018.11.043 .031 Ma L, Hernandez MO, Zhao Y, Mehta M, Tran B, Kelly M, Rae Z, Hernandez The Tabula Muris Consortium. 2020. A single-cell transcriptomic atlas char- JM, Davis JL, Martin SP, et al. 2019. Tumor cell biodiversity drives micro- environmental reprogramming in liver cancer. Cancer Cell 36: 418– acterizes ageing tissues in the mouse. Nature 583: 590–595. doi:10 430.e6. doi:10.1016/j.ccell.2019.08.007 .1038/s41586-020-2496-1 Mereu E, Lafzi A, Moutinho C, Ziegenhain C, McCarthy DJ, Álvarez-Varela Tauriello DVF, Palomo-Ponce S, Stork D, Berenguer-Llergo A, Badia- A, Batlle E, Sagar, Grün D, Lau JK, et al. 2020. Benchmarking single-cell Ramentol J, Iglesias M, Sevillano M, Ibiza S, Cañellas A, Hernando- RNA-sequencing protocols for cell atlas projects. Nat Biotechnol 38: 747– Momblona X, et al. 2018. TGFβ drives immune evasion in genetically re- 755. doi:10.1038/s41587-020-0469-4 constituted colon cancer metastasis. Nature 554: 538–543. doi:10.1038/ Merlos-Suárez A, Barriga FM, Jung P, Iglesias M, Céspedes MV, Rossell D, nature25492 Sevillano M, Hernando-Momblona X, da Silva-Diz V, Muñoz P, et al. Tirosh I, Izar B, Prakadan SM, Wadsworth MH, Treacy D, Trombetta JJ, 2011. The intestinal stem cell signature identifies colorectal cancer Rotem A, Rodman C, Lian C, Murphy G, et al. 2016a. Dissecting the stem cells and predicts disease relapse. Cell Stem Cell 8: 511–524. multicellular ecosystem of metastatic melanoma by single-cell RNA- doi:10.1016/j.stem.2011.02.020 seq. Science 352: 189–196. doi:10.1126/science.aad0501 Mimitou EP, Cheng A, Montalbano A, Hao S, Stoeckius M, Legut M, Roush Tirosh I, Venteicher AS, Hebert C, Escalante LE, Patel AP, Yizhak K, Fisher T, Herrera A, Papalexi E, Ouyang Z, et al. 2019. Expanding the CITE-seq JM, Rodman C, Mount C, Filbin MG, et al. 2016b. Single-cell RNA-seq tool-kit: detection of proteins, transcriptomes, clonotypes and CRISPR supports a developmental hierarchy in human oligodendroglioma. perturbations with multiplexing, in a single assay. Nat Methods 16: Nature 539: 309–313. doi:10.1038/nature20123 409–412. doi:10.1038/s41592-019-0392-0 Tumeh PC, Harview CL, Yearley JH, Shintaku IP, Taylor EJM, Robert L, Park J-E, Botting RA, Domínguez Conde C, Popescu D-M, Lavaert M, Kunz Chmielowski B, Spasic M, Henry G, Ciobanu V, et al. 2014. PD-1 block- DJ, Goh I, Stephenson E, Ragazzini R, Tuck E, et al. 2020. A cell atlas ade induces responses by inhibiting adaptive immune resistance. Nature of human thymic development defines T cell repertoire formation. 515: 568–571. doi:10.1038/nature13954 Science 367: eaay3224. doi:10.1126/science.aay3224 van der Wijst MGP, Brugge H, de Vries DH, Deelen P, Swertz MA, LifeLines Patel AP, Tirosh I, Trombetta JJ, Shalek AK, Gillespie SM, Wakimoto H, Cohort Study, BIOS Consortium, Franke L. 2018. Single-cell RNA se- Cahill DP, Nahed BV, Curry WT, Martuza RL, et al. 2014. Single-cell quencing identifies celltype-specific cis-eQTLs and co-expression RNA-seq highlights intratumoral heterogeneity in primary glioblasto- QTLs. Nat Genet 50: 493–497. doi:10.1038/s41588-018-0089-9 ma. Science 344: 1396–1401. doi:10.1126/science.1254257 van Dijk D, Sharma R, Nainys J, Yim K, Kathail P, Carr AJ, Burdziak C, Moon Peng J, Sun B-F, Chen C-Y, Zhou J-Y, Chen Y-S, Chen H, Liu L, Huang D, KR, Chaffer CL, Pattabiraman D, et al. 2018. Recovering gene interac- Jiang J, Cui G-S, et al. 2019. Single-cell RNA-seq highlights intra-tumoral tions from single-cell data using data diffusion. Cell 174: 716– heterogeneity and malignant progression in pancreatic ductal adeno- 729.e27. doi:10.1016/j.cell.2018.05.061 carcinoma. Cell Res 29: 725–738. doi:10.1038/s41422-019-0195-y Vieira Braga FA, Kar G, Berg M, Carpaij OA, Polanski K, Simon LM, Brouwer Popescu D-M, Botting RA, Stephenson E, Green K, Webb S, Jardine L, S, Gomes T, Hesse L, Jiang J, et al. 2019. A cellular census of human lungs Calderbank EF, Polanski K, Goh I, Efremova M, et al. 2019. Decoding hu- man fetal liver haematopoiesis. Nature 574: 365–371. doi:10.1038/ identifies novel cell states in health and in asthma. Nat Med 25: 1153– s41586-019-1652-y 1163. doi:10.1038/s41591-019-0468-5 Puram SV, Tirosh I, Parikh AS, Patel AP, Yizhak K, Gillespie S, Rodman C, Wu TD, Madireddi S, de Almeida PE, Banchereau R, Chen Y-JJ, Chitre AS, Luo CL, Mroz EA, Emerick KS, et al. 2017. Single-cell transcriptomic Chiang EY, Iftikhar H, O’Gorman WE, Au-Yeung A, et al. 2020. analysis of primary and metastatic tumor ecosystems in head and Peripheral T cell expansion predicts tumour infiltration and clinical re- neck cancer. Cell 171: 1611–1624.e24. doi:10.1016/j.cell.2017.10.044 sponse. Nature 579: 274–278. doi:10.1038/s41586-020-2056-8 Qian J, Olbrecht S, Boeckx B, Vos H, Laoui D, Etlioglu E, Wauters E, Pomella Yofe I, Dahan R, Amit I. 2020. Single-cell genomic approaches for develop- V, Verbandt S, Busschaert P, et al. 2020. A pan-cancer blueprint of the ing the next generation of immunotherapies. Nat Med 26: 171–177. heterogeneous tumor microenvironment revealed by single-cell profil- doi:10.1038/s41591-019-0736-4 ing. Cell Res 30: 745–762. doi:10.1038/s41422-020-0355-0 Yost KE, Satpathy AT, Wells DK, Qi Y, Wang C, Kageyama R, McNamara KL, Ramachandran P, Dobie R, Wilson-Kanamori JR, Dora EF, Henderson BEP, Granja JM, Sarin KY, Brown RA, et al. 2019. Clonal replacement of tu- Luu NT, Portman JR, Matchett KP, Brice M, Marwick JA, et al. 2019. mor-specific T cells following PD-1 blockade. Nat Med 25: 1251–1259. Resolving the fibrotic niche of human liver cirrhosis at single-cell level. doi:10.1038/s41591-019-0522-3 Nature 575: 512–518. doi:10.1038/s41586-019-1631-3 Zeisel A, Hochgerner H, Lönnerberg P, Johnsson A, Memic F, van der Zwan J, Regev A, Teichmann SA, Lander ES, Amit I, Benoist C, Birney E, Bodenmiller Häring M, Braun E, Borm LE, La Manno G, et al. 2018. Molecular archi- B, Campbell PJ, Carninci P, Clatworthy M, et al. 2017. Science Forum: tecture of the mouse nervous system. Cell 174: 999–1014.e22. doi:10 The Human Cell Atlas. eLife 6: e27041. doi:10.7554/eLife.27041 .1016/j.cell.2018.06.021 Rodriques SG, Stickels RR, Goeva A, Martin CA, Murray E, Vanderburg CR, Zhang Q, He Y, Luo N, Patel SJ, Han Y, Gao R, Modak M, Carotta S, Haslinger Welch J, Chen LM, Chen F, Macosko EZ. 2019. Slide-seq: a scalable tech- C, Kind D, et al. 2019. Landscape and dynamics of single immune cells nology for measuring genome-wide expression at high spatial resolu- in hepatocellular carcinoma. Cell 179: 829–845.e20. doi:10.1016/j.cell tion. Science 363: 1463–1467. doi:10.1126/science.aaw1219 .2019.10.003 Sade-Feldman M, Yizhak K, Bjorgaard SL, Ray JP, de Boer CG, Jenkins RW, Lieb DJ, Chen JH, Frederick DT, Barzily-Rokni M, et al. 2018. Defining T cell states associated with response to checkpoint immunotherapy in melanoma. Cell 175: 998–1013.e20. doi:10.1016/j.cell.2018.10.038 Received November 3, 2020; accepted in revised form August 17, 2021. 1926 Genome Research www.genome.org