Single-Cell Analysis of the Pan-Cancer Immune Microenvironment and scTIME Portal

The efficacy of immunotherapy is significantly influenced by the heterogeneity of tumor immune microenvironments (TIME; refs. 1–5), and a deep characterization of TIMEs will benefit the clinical applications. Generally, there are two types of TIMEs, “cold” and “hot.” The “cold” TIMEs are featured by a high density of tumor-associated macrophages (TAM), and the “hot” TIMEs are characterized by an increased infiltration of PDCD1 (PD-1)-expressing cytotoxic T lymphocytes (CTL; ref. 1). The PDCD1-high T cells are referred to dysfunctional/exhausted T cells, which have poor effector function (6–8). The intrinsic features of tumor cells trigger the distinct TIMEs, whereas the interactions between the immune cells in TIMEs also shape the complexity of TIMEs. It is well known that TAMs affect the T-cell responses in the tumor microenvironment (9). TAMs can inhibit effective adaptive immunity through multiple mechanisms, including promoting the function of regulatory T cells (Treg), metabolically starving the T cells, and triggering the T-cell expression of inhibitory immune checkpoints (10, 11). However, the common and tumor type–specific cross-talk among distinct immune cell subsets in TIMEs has not been fully studied in a pan-cancer view.

High-resolution characterization of TIMEs in various single tumor types was achieved by single-cell RNA sequencing (scRNA-seq) in an unprecedented manner (4, 6, 12–21), which largely extended our understanding of TIMEs. A high variation of TIME composition among patients is observed virtually in all studies (6, 12–19). For example, the fraction of dysfunctional T cells was observed to be highly variable in melanoma (6). Yet, it remains unclear whether the variation of the infiltration and accumulation of distinct immune cells in TIMEs are correlated across patients or different tumor types. Besides the abundance, the immune cell states in TIMEs are also heterogeneous across patients. For instance, infiltrating T-cell states are not discrete but manifest a continuous spectrum in activation and differentiation trajectory in breast cancer due to the substantial heterogeneity (16). Macrophages were shown to bear several subclusters, such as complement gene C1QC⁺ and secreted phosphoprotein gene SPP1⁺ subclusters in colon cancer, which are beyond the conventional M1 and M2 definition of macrophages (12–14). However, a single-cell pan-cancer survey of the heterogeneity and the tumor type–specific features of immune cell states in TIMEs has not been fully examined. This is particularly essential for understanding why the current immunotherapies only succeed in a subset of cancer types, such as in melanoma and lung cancer (22).

Here, we aimed to examine the common and specific immune features of TIMEs within and across tumor types to understand the complexity and heterogeneity of TIMEs better. We performed a comprehensive comparative analysis of the single-cell transcriptome of TIMEs across six human cancer types. With a total of 196,273 cells from 136 patients, we first built a pan-cancer TIME reference to coordinate the cell annotations across tumor types. We then characterized the heterogeneity of each cell type and the correlation of immune cells in TIMEs across patients and tumor types in terms of abundance, cell states, and cell communications.

With the thriving applications of single-cell technology, many single-cell databases have emerged (23–28), including Single Cell Portal (https://singlecell.broadinstitute.org/single_cell) and Single Cell Expression Atlas (https://www.ebi.ac.uk/gxa/sc/home). However, these databases mainly provide the visualization of cell clusters and gene expression levels. We built a single-cell transcriptome database and a Web portal of TIMEs, named “single-cell TIME” (scTIME), to facilitate the TIME-specific analysis at the single-cell level. With the single-cell decomposition of TIMEs, true immune signatures, rather than computational estimations, from bulk transcriptomes can be retrieved (29–31). All the human datasets in scTIME were annotated using the pan-cancer reference generated from this study to make the comparative analysis across datasets possible. We envision the scTIME portal will benefit the cancer researchers and clinicians, and together with our pan-cancer analysis, will deepen our understanding of the complexity of TIMEs and immunotherapeutic outcomes across patients or tumor types at a single-cell level.

This study does not involve any patient or animal samples. Data used in this study were publically available and downloaded from the Gene Expression Omnibus and Genome Sequence Archive as shown in Supplementary Table S1. The datasets were selected based on two criteria. First, the datasets were selected from articles that had been peer-reviewed. Second, to explore the full features and relationship across different cell types of TIME, we mainly chose the datasets that contained both lymphoid and myeloid lineages and more than 500 cells of each lineage. The gene quantification matrix of each dataset was downloaded. All genes were transferred to Ensembl v100 using R package biomaRt. Only cells with expression of CD45 or CD3E (counts ≥ 1) and mitochondria gene counts lower than 30% were kept for the downstream analysis. After filtering, 7.6% of cells were CD45⁻CD3E⁺.

Due to the inherent feature of the sequencing technologies and platforms, doublets are mostly excluded in the experimental procedures of SMART-Seq2 and massively parallel scRNA-seq (MARS-seq; ref. 6). The doublet rate of the inDrop method used in the Breast dataset was approximately 0.59% (16). Thus, we only performed the doublet identification for the Lung dataset (inDrop platform) and 10X platform datasets of Colon, Liver2019, and pancreatic ductal adenocarcinoma (PDAC) using Scrublet (https://github.com/AllonKleinLab/scrublet; ref. 32). The expected doublet rate was set to 6%, which is suggested by 10X Chromium (10X Genomics). We identified 745 (0.593%) doublets from the Colon, Liver2019, Lung, and PDAC datasets (Supplementary Fig. S1A). Due to the low doublet count, we did not filter the doublets identified for the downstream analysis.

Datasets were integrated (batch effect corrected) at normalized count level across patients, tissues, and platforms with the functions of FindIntegrationAnchors and IntegrateData from Seurat v3 (https://satijalab.org/seurat/; ref. 33). The integrated data were then clustered by Seurat using the Louvain algorithm with resolution 2, and the myeloid cells were extracted and reclustered with resolution 1. The cell type annotations were determined based on the marker genes of each cluster identified by Seurat and the original cell annotation from the associated publications. Endothelium and fibroblast cells were excluded, which weakly expressed CD45 or CD3E (median counts of 2.5 and 0). Cluster dendrograms were built with the top 60 principal components (PC) of 2,000 feature genes (i.e., most highly variable genes) of the integrated data. The top 60 PCs were selected based on both the variance elbow plot and the significance of the PCs estimated by the JackStraw procedure from Seurat. 23.5% variance was explained by the top 60 PCs. A cell type was defined as enriched in a certain tissue if the ratio of observed to expected fraction of the cell type was > 1.5 (χ² test, P < 0.001). M1 and M2 signature genes were downloaded from Azizi and colleagues (16).

We employed SCENIC (https://github.com/aertslab/SCENIC) to analyze the activity of transcription factors (TF; ref. 34). In brief, SCENIC groups the target genes regulated by a TF as a regulon and calculates the activity of the regulon as the relative rank-sum of the expression of these targets. Myeloid cells and T/natural killer (NK) cells were subsampled to 10,000 each, and the regulon analysis was performed separately. To preserve rare cell clusters, cells from tumors were subsampled using a geometric sketching method (35). Finally, the binary activity of regulons was estimated by SCENIC (34). The cell type–enriched regulons were the regulons activated in >50% cells of a given cell type, regardless of tumor type, but not activated in >80% cells of any single tumor type. The regulons enriched in all cell types of myeloid or T/NK cells were excluded.

For cellular composition analysis, only the datasets with CD45⁺ sorted cells or red blood cell (RBC)–depleted cells (Supplementary Fig. S1B) and the cells expressing CD45 were used. For the correlation analysis of cell fraction, treatment-naïve patients with more than five cells in at least three cell types were kept. Cell types with more than 5 cells in at least 5 patients with a cell fraction ≥0.5% in at least 10 patients were selected. For the variation analysis of the cell fraction, all cell types were included. Cell fractions <0.5% were assigned to 0.5% to avoid infinity values after log transformation. The median and variance of cell fractions were calculated using the log-transformed values.

The integrated data by Seurat introduced dependencies between data points and violated the assumptions of the statistical tests used for differential expression (https://github.com/satijalab/seurat/discussions/4000). To correct the batch effect for the analysis of transcriptome heterogeneity and differential expression, we generated pseudo bulk samples and used the function removeBatchEffect from R package edgeR (36) to remove batch effects. The original transcriptome of the cells from the same patient, the same cell type, and the same platform were averaged into one pseudo bulk sample. A total of 3,385 pseudo samples from the tumor samples were generated. We treated the global difference of all cell types across patients and platforms as batch effects and kept the cell type–specific difference among samples. Genes that were commonly detected in all eight datasets were used. Three of the 136 patients we compiled had samples profiled on two different sequencing platforms, thus, after the batch effect correction, samples from different sequencing platforms, but from the same patient and cell type, were averaged into one sample. Samples that were averaged from less than three cells were filtered out, and finally, 1,940 treatment-naïve samples were selected for the downstream analysis.

The variance of the transcriptome was estimated in the Uniform Manifold Approximation and Projection (UMAP) reduction space (UMAP reduction performed using Seurat) of the pseudo samples. Outliers based on the boxplot statistics of each cell type were filtered out. Outliers were samples with a value greater than Q3 + 1.5 × IQR or less than Q1 – 1.5 × IQR. Q3 is 75th percentile, Q1 is 25th percentile, and IQR is the interquartile range. R function boxplot.stats were used to calculate the threshold of the outliers. The highly variable genes of CD8-PDCD1 were estimated by method “vst” in Seurat v3. Enrichment analysis of the top 500 highly variable genes was performed by Metascape (http://metascape.org; ref. 37).

In each cell type, differentially expressed genes (DEG) between one cancer type and all other cancer types were detected with the pseudo bulk samples using R package limma (38). DEGs were defined as fold change >1.5 and adjusted P < 0.05. Enrichment analysis of the top 500 DEGs was performed by R package clusterProfiler (39). Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enriched with gene counts ≥8 and adjusted P value < 0.05 were selected, and the pathways related to immunity are shown in figures.

The datasets with FACS-purified CD45⁺ cells or RBC-depleted cells were used. Treatment-naïve patients with more than five cells in at least three cell types were included. Cell types with more than 5 cells in at least 5 patients and a cell fraction ≥0.5% in at least 30 patients were analyzed. In each cell type, pseudo bulk samples averaged from more than or equal to three cells were selected. PCs of the top 2,000 highly variable genes of the pseudo bulk transcriptomes were defined as gene signatures in each cell type. The highly variable genes were estimated by method “vst” in Seurat v3. Pairwise Spearman correlations were calculated between all PCs of the transcriptome of one cell type and the fraction of another cell type. Outliers of the signatures based on the boxplot statistics were filtered out. Enrichment analysis of the genes weighted top 30 in the PC was performed by Metascape.

Ligand–receptor (LR) network analysis was performed using CellPhoneDB v2 (https://www.cellphonedb.org/; ref. 40) with the parameters: “–iterations = 1000 –subsampling –subsampling-log false –subsampling-num-cells 1000 –threshold 0.01.” Predicted LR pairs with P < 0.01 were selected. To be comparable across datasets, we normalized the LR mean strength. The mean was normalized to the total strength predicted in each dataset with P < 0.01 and then multiplexed by 10⁴, that is mean strength per 10⁴ total mean strength. The LR pairs with neither partner annotated as receptors were filtered out. The normalized counts of ligand or receptor were the counts per 100 total counts of ligand and receptor that formed significant interaction pairs in the dataset.

The VueJS v2.6.0 (https://vuejs.org/) and ViewUI v4.0.0 (https://www.iviewui.com/) were used to develop the server front end, whereas the back end was built in Java using the SpringBoot Web framework v2.1.13. The server is hosted on a Linux Centos v7.0.1406 server running Apache Tomcat v9.0.31. The files, including raw data matrix, preprocessed data, and intermediate results, were stored in the local Ext4 file system, whereas the metadata and the analysis results were stored in the MySQL v.8.0.18 database. This enables more efficient querying, visualizing, and archiving of the datasets and analysis results. The interactive visualization diagrams are implemented with the D3.js (https://d3js.org/), Echarts.js (https://echarts.apache.org/), and CanvasXpress.js (http://canvasxpress.org/). Except for the interactive visualization of the Web server, all other steps of data analysis are performed by the corresponding R functions.

We first searched PubMed for literature with scRNA-seq related to the cancer immune microenvironment. We filtered and collected 49 available datasets with immune cell numbers more than 500 or with original cell labels. We then curated the metainformation of each dataset, including species, cancer type, technology, tissue source, flow cytometry gates, perturbations, and so on. We downloaded the quantification matrix of gene expression. Each dataset was lifted over to Ensembl v100 using R package biomaRt, normalized, and integrated across samples using function FindIntegrationAnchors and IntegrateData from Seurat v3. For human samples, the immune cells in each dataset were reannotated with the cell type reference generated in this study using function FindTransferAnchors and TransferData of Seurat v3. Curated metainformation and processing steps of each dataset were elaborated in the database.

In the LR interaction network module, we use CellPhoneDB v2 to predict the communications between cell types, using the parameters: “–iterations = 1000 –subsampling –subsampling-log false –subsampling-num-cells 1000 –threshold 0.01.” The mean strength was normalized to the total strength predicted in each dataset at a user-defined P value cutoff and then multiplied by 10⁴.

scTIME is a database and a Web portal available at http://scTIME.sklehabc.com.

To pave the way for the pan-cancer comparative analysis of TIMEs at a single-cell level, we first built a pan-cancer immune cell reference by integrating eight single-cell transcriptome datasets across six types of cancers (Supplementary Table S1). The cancer types include breast carcinomas (16), non–small cell lung cancer (12), colorectal cancer (14), hepatocellular carcinoma (liver cancer; refs. 13, 17), melanoma (6, 41), and PDAC (19). A total of 196,273 cells from 136 patients passed the quality controls (Supplementary Table S1). We integrated the datasets to correct the batch effects across patients or sequencing technologies (Supplementary Fig. S1C and S1D) for cell clustering using Seurat and obtained 45 consistent immune clusters among all cancer types (Fig. 1A). We renamed the clusters with the marker genes identified from the pan-cancer dataset (Fig. 1B and C; Supplementary Table S2).

We obtained 9 CD4⁺ T-cell clusters, including naïve T cells (CCR7⁺, including TCF7⁺, IL7R⁺, and FOS⁺ subsets), memory T cells (IL7R⁺, including an MAL⁺ subset), T follicular helper cells (Tfh; CXCL13⁺), and Tregs (Fig. 1A and B). Tregs consisted of CTLA4⁺ and LEF1⁺ subclusters, representing effector and resting Tregs, respectively (42). CD8⁺ T cells were partitioned into nine clusters, including naïve T cells (CCR7⁺), mucosal-associated invariant T cells (MAIT; SLC4A10⁺), memory T cells (GZMK⁺, including a FOS⁺ subset), predysfunctional T cells (ZNF683⁺; refs. 15, 17), dysfunctional T cells (PDCD1⁺), and cytotoxic T cells (GNLY⁺; Fig. 1A and B; Supplementary Fig. S1E). CTLA4⁺ Tregs, CD8-HSPA1A, and CD8-PDCD1 were identified to be enriched in tumors (see Materials and Methods). CD8-ZNF683, MAIT cells, and CD8-GZMK-FOS were prevalent in tumor adjacent normal tissue or normal donors, suggesting that they are likely tissue-resident immune cells (Fig. 1D; Supplementary Fig. S2A). Two naïve T-cell clusters (CD4-CCR7-TCF7 and CD8-CCR7), CD4-IL7R-MAL, CD4-LEF1-Treg, and cytotoxic CD8-GNLY were prevalent in peripheral blood, whereas other naïve T cells (CD4-CCR7-IL7R and CD4-CCR7-FOS) and memory T cells (CD4-IL7R and CD8-GZMK) were prevalent in lymph nodes (Fig. 1D; Supplementary Fig. S2A). We identified an innate lymphoid cell cluster and three NK-cell clusters expressing KLRD1 (Fig. 1A and B), including the classic FCGR3A⁺ NK cells, NK-KLRC1, and NK-ITGA1. The dendrogram of the lymphoid clusters shows the transcriptomic similarity across the clusters (Fig. 1E).

Eight clusters were identified as TAMs based on the expression of CD68 and CD163 (Fig. 1A and C). Two C1QC⁺ TAM clusters (C3⁺ and PLTP⁺, respectively) expressed a high level of complement C1Q genes and CSF1R (Fig. 1C), indicating its phagocytosis function and high sensitivity to anti-CSF1R treatment (14). Both SPP1⁺ TAM clusters expressed a higher SPP1 and MARCO compared with C1QC⁺ TAMs. Macrophage-SPP1-CLEC5A expressed higher angiogenesis genes VEGFA and EREG and inflammation-related gene CLEC5A, whereas Macrophage-SPP1-ACP5 expressed higher TAM markers APOE, MRC1, and TREM2 (Fig. 1C). All the C1QC⁺ and SPP1⁺ TAMs were biased toward an M2 signature (Supplementary Fig. S2B; Fig. 1D). Macrophage-CXCL10 had high expression of the chemokines CXCL10 and CXCL9, which regulate T-cell infiltration (Fig. 1C; Supplementary Table S2; refs. 43, 44). Macrophage-IL1B had high expression of myeloid-derived suppressor cell signature genes FCN1, VCAN, IL1B, THBS1, NLRP3, and VEGFA and more resembled monocytes and neutrophils than other macrophage clusters (Fig. 1C and F; refs. 13, 45). Macrophage-CXCL10 and Macrophage-IL1B were biased toward an M1 signature (Supplementary Fig. S2B). Macrophage-MARCO had high expression of the TAM marker gene MARCO, a scavenger receptor (Fig. 1C; ref. 46). Macrophage-AREG expressed AREG and lower level of other macrophage marker genes than M2- and M1-biased TAMs, indicating an intermediate state (Fig. 1C). Four dendritic cell (DC) types were identified, including conventional DCs (cDC1 and cDC2), DC-LAMP3 (13), and plasmacytoid DCs (pDC; Fig. 1A). Tumors were enriched for M2-biased TAMs, as well as Macrophage-CXCL10, DC-LAMP3, and pDCs (Fig. 1D). Normal tissue and tumors showed comparable enrichment for Macrophage-AREG, Macrophage-IL1B, and cDC2. Macrophage-MARCO and cDC1 were prevalent in tumor adjacent normal tissue or normal donors and rare in the blood or lymph nodes, suggesting that they are likely tissue-resident immune cells (Fig. 1D).

We also used SCENIC (ref. 34; see Materials and Methods) to explore the cell type–enriched transcriptional regulatory network in the pan-cancer dataset. The TFs PPARG, MAF, MAFB, NR1H3, and PBX3 showed similar regulon activities across the four M2-biased TAM clusters, whereas the regulon activities of MITF, BHLHE41, CEBPA, and ZNF627 were higher in Macrophage-SPP1-ACP5 than in Macrophage-SPP1-CLEC5A (Fig. 1G). ETV5 was more active in C1QC⁺ TAMs (Fig. 1G). Macrophage-CXCL10 was enriched with higher regulon activity of ETV7 and HIVEP3 (Fig. 1G). In addition to the previously reported BCL3 and RXRA (13), we found the regulon activities of BACH1, CREB5, NFE2L2, NFIL3, and TGIF1 were also high in Macrophage-IL1B (Fig. 1G). Differential activation of these TFs may underlie the distinct functional regulation of these TAMs in the TIME. The cell type–enriched activated TFs were also observed in T cells and NK cells, such as LEF1, TCF7, FOXP3, and EOMES, as reported in previous studies of single cancer types (Supplementary Fig. S2C; refs. 6, 14, 15).

To map conserved and tumor type–specific features of the TIMEs, we analyzed the TIME composition of treatment-naïve patients in each tumor type. For consistency across different studies, we chose datasets generated from FACS-purified CD45⁺ cells or RBC-depleted cells (Breast, Liver2019, Lung, Melanoma2018, and PDAC) and the cells with CD45 expression for the composition analysis of TIMEs. The composition of the TIMEs was heterogeneous across patients and tumor types (Fig. 2A; Supplementary Fig. S3), as observed in the previous studies (6, 12–19). We found that the fractions of the dysfunctional CD8⁺ T cells, TAMs, and cDC2 were highly variable across the patients in TIMEs (Fig. 2A), similar to that observed in melanoma (6). High variation in the abundance of these cells may underlie disease severity and patient survival, as a high marker gene expression of dysfunctional CD8⁺ T cells and TAMs associates with poor prognosis of patients with liver cancer (13, 15); however, cDC2 marker gene expression manifested a positive association with the survival of patients with lung cancer (12). In contrast, we found that the fraction of CTLA4⁺ Tregs, which were also enriched in tumors (Fig. 1D), was relatively stable (Fig. 2A). This indicated that CTLA4⁺ Tregs might be commonly involved in T-cell inhibition across patients and tumor types.

Given the cross-regulation of different immune cell types (1, 10, 47, 48), we explored the pairwise correlation of the cell fractions of immune cell types in TIMEs across patients. The fractions of M2-biased TAMs were inversely correlated with the fractions of Tfh and CD8-PDCD1, both of which had high expression of PDCD1 (Figs. 1B and 2B; Supplementary Fig. S4A). This is in agreement with the observations that in “cold” tumors, high TAM abundance decreases T-cell infiltration, and “hot” tumors with low TAM abundance and high T-cell infiltration are featured by PDCD1-high CTLs (1). However, there was no significant correlation between the fractions of TAMs and cytotoxic CD8-GNLY. The fraction of CD8-PDCD1 was positively correlated with Tfh (Fig. 2B), as observed in melanoma (6), but the correlation with naïve T cells was not significant.

Although the composition of TIMEs varies across patients, it was more similar across patients within the same tumor type compared with patients from different tumor types (Fig. 2C; Supplementary Figs. S3A–S3C and S4A–S4C). Thus, the high variations of TIME composition were partially due to the difference between tumor types (Fig. 2D; Supplementary Fig. S4B). For example, the fraction of CD8-PDCD1 exhibited low variance across patients within each tumor type (Fig. 2D; Supplementary Fig. S4B), but the mean fractions of CD8-PDCD1 varied substantially across different tumor types (Fig. 2D). By comparing TIMEs across tumor types, we found the majority of the cell types were shared across tumor types, although the abundance of many cell types showed different biases in different tumor types (Fig. 2E; Supplementary Fig. S4C). Compared with the other tumor types, lung cancer and PDAC presented a higher fraction of macrophages (24.1% and 26.8%, respectively) and lower T-cell infiltration (Supplementary Fig. S3D). Notably, the abundance of macrophages in normal pancreatic tissue was also high (58.7%; Supplementary Fig. S3D). Further examination showed that lung cancer and PDAC were enriched for M2-biased TAMs and a subset of naïve CD4⁺ T cells (CD4-CCR7-TCF7 and CD4-CCR7-IL7R; Fig. 2E; Supplementary Fig. S3E and S3F). Melanoma was enriched for the CD8-PDCD1 cells, whereas breast and liver cancer were enriched for cytotoxic CD8-GNLY, CD8-ZNF683, and NK cells (Fig. 2E; Supplementary Fig. S3F). Thus, our findings suggest that high CD8-PDCD1 abundance may underlie why patients with melanoma generally respond better to immune checkpoint therapy than patients with non–small cell lung cancer or PDAC (49). These composition biases of the TIMEs were also observed at a single patient level (Supplementary Fig. S3A–S3C). However, the difference in the immune cell abundance across tumor grades or metastasis stages was not significant for most cell types (Supplementary Fig. S4D).

Not merely the difference in immune cell abundance, the transcriptional variance of the immune cells also shapes the patient-specific TIMEs. We next examined the transcriptional variation of each immune cell type across patients. Because the integrated data introduced dependencies between datapoints and violated the assumptions of the statistical tests used for differential expression, we created pseudo bulk samples to perform the batch effect correction as bulk RNA-seq for the analysis of transcriptome heterogeneity and differential gene expression (see Materials and Methods). We averaged the gene expression across the cells of each cell type in each patient to create pseudo bulk samples (Supplementary Fig. S5A and S5B). To examine the transcriptome variance of each cell type across patients or tumor types, the global difference across patients was removed as batch effects.

Considering the population size of each cell type, CD8-PDCD1 showed a high transcriptome variance across patients (Fig. 3A; Supplementary Fig. S5C), whereas CTLA4⁺ Tregs and Tfh were relatively stable (Fig. 3A). This is in line with a dynamically differentiating and active proliferative state of the dysfunctional T cells in TIMEs (6, 50). The highly variable genes of CD8-PDCD1 cells were enriched for the functions related to cytokine interaction and interaction between lymphoid and nonlymphoid cells (Fig. 3B; Supplementary Fig. S5D). Compared with CD8-PDCD1, CTLA4⁺ Tregs exhibited relatively stable features in both cell fraction and cell state (Figs. 2A and 3A), hinting that Tregs may be regulated by conserved mechanisms across patients and tumor types.

Although it is known that immune cells form a complex regulatory network in TIMEs, the precise relationships between the subtypes of immune cells are not clear. To take advantage of the refined annotation of immune cell types in the TIME, we in silico screened the transcriptional signatures of one immune cell type that correlated to the abundance of other immune cell types. The PCs of the top 2,000 highly variable genes in each cell type were used as the transcriptional signatures. Focusing on the correlation between myeloid cells and T cells with coefficient |R|≥0.6, we found the signatures of the M2-biased TAMs, Macrophage-AREG, classical monocyte (Monocyte-VCAN), and cDC2 were correlated with the fraction of CD8-PDCD1 and memory T cells (CD8-GZMK and CD4-IL7R; Supplementary Fig. S6A; Supplementary Table S3). In particular, a signature of the PLTP⁺ subset of C1QC⁺ macrophages (Macrophage-C1QC-PLTP), comprising TNFSF13, CCL2, CXCL12, SLC40A1, etc., positively correlated with the fraction of CD8-PDCD1 (Fig. 3C and D). This signature of Macrophage-C1QC-PLTP was enriched for functions related to chemokine-mediated signaling and regulation of complement, exocytosis, extracellular matrix (ECM), and cell adhesion (Fig. 3C and D), suggesting these biological events in Macrophage-C1QC-PLTP may be involved in the regulation of the dysfunctional T-cell abundance in TIMEs. SLC40A1 is reported to be associated with the prognosis of liver cancer (13). Our findings suggest that the prognostic association of SLC40A1 may be established through regulating the fraction of CD8-PDCD1 by TAMs. A signature of cDC2, including genes involved in leukocyte activation and necroptosis, negatively correlated with the fraction of CD8-PDCD1 (Fig. 3E and F). On the other hand, the gene signatures of CD8-PDCD1, Tregs, and Tfh were correlated with the fractions of myeloid cells, especially cDC2 (Supplementary Fig. S6A–S6D), suggesting a reciprocal regulation between the myeloid and T cells.

Next, we asked whether there were transcriptomic features of TIMEs specific to certain tumor types. By comparing each tumor type with the others, we recognized the DEGs of each cell type across tumor types (Supplementary Table S4). KEGG pathway analysis of significantly upregulated genes in the M2-biased TAMs and Macrophage-AREG in lung cancer revealed their participation in cytokine interaction, chemokines, and TNF signaling pathways (Fig. 3G and H; Supplementary Fig. S6E; Supplementary Table S5), especially for proinflammatory pathways (Supplementary Fig. S6F). Thus, although both lung cancer and PDAC exhibited high TAM abundance, M2-biased TAMs in lung cancer might have elevated proinflammatory function (Fig. 2E; Supplementary Fig. S6F). Genes involved in the PI3K–Akt signaling pathway and focal adhesion were significantly upregulated in cDC2 in PDAC (Fig. 3G). Genes significantly downregulated in M2-biased TAMs in PDAC were enriched for cytokine interaction and chemokine signaling pathways (Supplementary Fig. S6G). These tumor type–specific expression features of certain immune cells may help design tumor type–oriented immunotherapies.

Cytokine/chemokine signaling plays a key role in regulating the recruitment and function of immune cells in TIMEs (9). However, a systematic quantification and comparison of the immune cell communication network across tumor types is lacking. We investigated immune cell communication mediated through LR interactions using CellPhoneDB (40). To compare across cancer types and scRNA-seq platforms, the mean strength of LR interactions was normalized within each dataset (see Materials and Methods). In general, the predicted pairwise interaction strength was strongest among macrophages and weakest among the lymphoid cells (Fig. 4A). TAMs were the core of the interaction network, which has been similarly observed in colon cancer (14). Because the close relationship between the macrophages and T cells was noted above, we then surveyed the interactions between macrophages and T-cell subsets. Macrophages were predicted to mainly communicate with CTLA4⁺ Tregs, Tfh, CD4-IL7R, CD8-GZMK, CD8-PDCD1, and CD8-ZNF683 T cells (Fig. 4A; Supplementary Fig. S7A). Among the macrophages, M2-biased TAMs, Macrophage-AREG, and Macrophage-CXCL10 communicated most strongly with T cells and NK cells (Fig. 4A; Supplementary Fig. S7A). Macrophage-CXCL10 and Macrophage-C1QC-PLTP were predicted to emanate stronger and more diverse signals to T cells and NK cells compared with other macrophages, whereas CTLA4⁺ Tregs were predicted to receive stronger and more diverse signals from macrophages rather than other T cells or NK cells (Fig. 4A; Supplementary Fig. S7A). These findings indicate that tumor-enriched immune cell types, such as CTLA4⁺ Tregs, CD8-PDCD1, and M2-biased TAMs, are the dominant cells in the interaction network between macrophages and T cells.

We next compared the immune cell interactions across tumor types. The split view of LR interactions in each tumor type displayed tumor-specific patterns (Supplementary Fig. S7B). The normalized counts of ligands and receptors were higher in macrophages from lung cancer, breast cancer, and PDAC (Fig. 4B), which mostly contributed to the myeloid–myeloid interactions (Fig. 4C, purple). The predicted LR interaction pairs between myeloid and lymphoid cells were similar across cancer types (Fig. 4C). We then further examined the communication between myeloid and lymphoid cells. Among T cells, we found macrophages exhibited strong interactions with CTLA4⁺ Tregs in all cancer types (Fig. 4D), as previously observed in liver cancer (20). This could underlie the stable features of the abundance and transcriptional state of CTLA4⁺ Tregs across tumor types (Figs. 2A and 3A). In contrast, among the other T-cell subsets, macrophages preferred different interaction partners in different tumor types (Fig. 4D). For instance, macrophages in colon cancer were predicted to interact more strongly with CTLA4⁺ Tregs and memory T cells (CD4-IL7R and CD8-GZMK; Fig. 4D, top three blue dots), whereas in liver cancer, macrophages were predicted to interact more strongly with CTLA4⁺ Tregs, CD8-PDCD1, and Tfh (Fig. 4D, top three green dots) compared with other T cells. On the other hand, among macrophages, T cells interacted most strongly with Macrophage-C1QC-PLTP (Supplementary Fig. S8A). Similarly, CTLA4⁺ Tregs and cDC2 were mostly involved in the communication between DC and T cells (Supplementary Fig. S8B and S8C). These findings indicate that the variations of cell interaction strength might shape different TIMEs in different tumor types.

We further focused on the communications among the tumor-enriched T cells and myeloid cells. Among 288 pairs of T-cell and macrophage subsets, the LR interaction strength of 17 pairs (5.9%) significantly correlated with their cell fractions (P < 0.05; Supplementary Table S6). The fraction of CD8-PDCD1 in T cells positively correlated with its predicted interaction strength with Macrophage-C1QC-PLTP, but not with other TAM subsets (Fig. 4E; Supplementary Fig. S8D; Supplementary Table S6). The LR pairs, including CD274 (PD-L1)-PDCD1, PDCD1LG2 (PD-L2)-PDCD1, CXCL9/10/12-CXCR3, TNFSF9-TNFRSF9, and TNFSF13-TNFRSF14, were generally enriched in colon cancer, liver cancer, and melanoma (Fig. 4F), in which the abundance of CD8-PDCD1 were higher (Fig. 4E). The function of these LRs was enriched for negative regulation of immune system process, TNF binding, and positive regulation of cell death (Fig. 4G). These findings again suggest that the increased cross-talk between the PLTP⁺C1QC⁺ TAMs and dysfunctional CD8-PDCD1 might promote the accumulation of the dysfunctional T cells. We also found the fraction of the predysfunctional T cell (CD8-ZNF683) and Tfh positively correlated with their LR interaction strength with Macrophage-SPP1-CLEC5A and Macrophage-MARCO, respectively (Supplementary Fig. S8E and S8F). These results suggest that the differential interaction strength with different TAM subsets may be responsible for the accumulation of different T-cell subsets in TIMEs.

Finally, to make full use of these single-cell transcriptome datasets of the TIMEs and facilitate the experimental researchers and clinicians to explore the features of TIMEs, we built the scTIME Portal (http://scTIME.sklehabc.com). The scTIME Web portal is a database for exploring and analyzing TIMEs at a single-cell level with the TIME-specific analysis modules (Fig. 5A; Supplementary Fig. S9A and S9B). scTIME allows the users to perform similar analysis that we performed above for each dataset on the portal. We collected scRNA-seq datasets of TIMEs from the Gene Expression Omnibus, Single Cell Portal, Single Cell Expression Atlas, and Genome Sequence Archive. The portal contains 49 datasets from 46 publications and includes 39 cancer types from human and mouse at present (Supplementary Table S7). We performed quality control, batch correction, normalization, and manually adjusted the metainformation for each dataset (Fig. 5A). Users can filter or search datasets, and then explore or analyze the dataset within the portal using the TIME-specific analysis modules (Supplementary Fig. S9B).

Five TIME-specific analysis modules were developed beside the UMAP exhibition of cell information and gene expression (Fig. 5A and B; Supplementary Fig. S9C). All the modules are interactive and user-friendly. The first and second “Composition” modules show immune cell composition of the TIME in each patient and the pairwise correlation of the cell fraction across cell types, respectively (Supplementary Fig. S9D and S9E). The third “Signature” module presents the relative expression of genes in selected signatures, such as the TNF signaling pathway and PD-1 checkpoint pathway (Fig. 5C). It can be used to compare the signature expression either across patients in a selected cell type or across cell types. It provides the information on the signature variance across patients or the signature enrichment across cell types. The fourth “Correlation” module allows users to explore the correlation between the gene expression in a cell type and the cell fraction of another cell type across patients (Fig. 5D). This module is designed to explore the genes that may regulate the composition of TIMEs. For example, users can explore whether CSF1R expression in the PLTP⁺C1QC⁺ TAMs is correlated with the fraction of dysfunctional CD8⁺ T cells in each dataset. The fifth “LR network” module exhibits the predicted LR network of selected cell types (Fig. 5E). In this module, users can select a subset of cell types, select network layout (Force-Directed and Circle), or filter edge strength to customize and explore the LR network. By switching the analysis type to “LR pair strength,” users can also view the strength of LR interaction pairs between two cell types (Supplementary Fig. S9F).

Another feature of scTIME is the unified cell annotation across datasets. Many downloaded datasets do not contain cell type annotations, which hinders the researchers from exploring the published data. We provide two versions of cell type annotations for all human datasets. One is the original annotation and the other is the unified annotation using the reference generated in this study (Fig. 5A). Users can switch the annotations easily with a slider button (Fig. 5B). In addition, we retain all types of cells, including malignant cells and fibroblasts, in our database. Users can also perform the above analysis beyond immune cells.

All the functions of these modules can be customized using specific filter parameters, including tissue type, patient phenotype, flow cytometry gating, treatment, and so on. Users can select a specific group of cells and samples of interest to analyze. For instance, users can select flow cytometry–sorted CD3⁺ cells from tumor and blood samples of treatment-naïve patients for subsequent analysis and comparison. The label transfer, cell fraction, correlation analysis, and LR network are precalculated for each dataset, allowing fast access to the results in these modules.

Extensive studies of TIMEs using scRNA-seq tend to be focused on individual cancer type, and some include a limited comparison with one or two of other tumor types (6, 12–19). Here, we systematically compared TIMEs from six cancer types in terms of cellular composition, cell states, and cell communications. The pan-cancer analysis enabled us to delicately identify cell subsets with increased cell numbers and examine the correlation of cell abundance and gene signatures with an increased number of patients.

It was observed that both cellular composition and cell states of TIMEs substantially vary across patients, which are influential factors to the different outcomes of immune therapies (4, 6, 12–19). Our pan-cancer analysis revealed that the fraction and cell state of the dysfunctional CD8⁺ T cells are highly variable, whereas those of CTLA4⁺ Tregs were relatively stable. This difference suggests that dysfunctional CD8⁺ T cells may have more critical roles in the differential response of distinct tumor types to immune therapies compared with CTLA4⁺ Tregs, even though both of them are enriched in tumors. Moreover, our pan-cancer analysis of cell communications in TIMEs demonstrated that the predicated strength of LR interactions varied across tumor types, especially for the interaction strength between TAMs and T cells. Among T cells, TAMs interacted most strongly with Tregs in all tumor types. Together with the relative stability of Tregs and high variation of CD8-PDCD1 in both abundance and transcriptional states, these findings suggest that targeting Tregs may have a broader impact on patient prognosis than targeting CD8⁺ T cells across distinct tumor types (51).

TAMs can trigger inhibitory immune checkpoints in T cells (10). A subset of TAMs (C1QC⁺ TAMs) was found to have high expression of genes related to PD-1 signaling and to interact with T-cell subsets in colon cancer (14). Through our pan-cancer analysis, we found the fraction of dysfunctional CD8⁺ T cells may be regulated by a specific subset of TAMs that express higher PLTP and C1QC through cytokine-mediated signaling. We further confirmed this association through LR interaction analysis. Meanwhile, our pan-cancer analysis revealed an inverse correlation between the fraction of C1QC⁺ TAMs and dysfunctional CD8⁺ T cells. We presumed that TIMEs with fewer TAMs have more T-cell infiltration, and that the Macrophage-C1QC-PLTP may trigger the conversion of infiltrating T cells into a dysfunctional state. High infiltration of TAMs reduces T-cell infiltration and thus decreases the net abundance of dysfunctional T cells. Although Macrophage-C1QC-PLTP could trigger inhibitory immune checkpoint expression on T cells, the low amount of infiltrating T cells may limit the increase of dysfunctional T cells. Based on the preferential depletion of the C1QC⁺ TAM population after anti-CSF1R treatment (14), we speculate that anti-CSF1R treatment will promote T-cell infiltration (52) and expand the dysfunctional T-cell population, which will provide more potential targets for immune checkpoint blockade therapies and improve their efficacy.

In addition to the immune cells, the molecular features of tumor cells, cancer-associated fibroblasts (CAF), and the ECM orchestrate a highly variable TIME (9, 48). The perturbations of these factors provide a possibility to unveil the robust correlation between immune cells in TIMEs. In this study, we focused on the tumor-enriched macrophages and T cells, whereas CAFs and tumor cells were not considered. CAFs and tumor cells are included in our database, and their roles in TIMEs can be explored using the scTIME Portal.

scTIME is a single-cell database and Web portal with TIME-specific analysis modules. Compared with TISCH (28) and other general single-cell databases (23, 24, 26, 27), scTIME provides TIME-specific analysis modules and a cell communication module. Similar to the bulk TIME analysis tools, such as TIMER (29) and CIBERSORT (30), we provide analysis for immune cell composition and correlation, but with precise cell type classification at the single-cell level. Cell communication and cell type–specific gene signature analysis are also included, which are specific for single-cell transcriptome analysis. Compared with the bulk TIME analysis database, a drawback of scTIME is the lack of survival analysis due to the limited number of patients. However, with the increasing amount of single-cell studies of TIMEs and the availability of patient survival information, this functional module can be added to scTIME in the future. We also unified the cell annotation, which allows comparison across datasets. Another single-cell database, Sfaira, also provides unified annotations with different levels of coarseness using pretrained models, although the immune cell annotations are still relatively rough (53). Sfaira mainly provides automated basic analysis of dimensionality reduction, clustering, and annotation using pretrained models, whereas scTIME provides TIME-specific analysis modules for deep exploration of the TIME. Svensson and colleagues provided curated descriptions of published single-cell datasets using Google Sheets, which is a useful source for data searching (54). Together with other publically available, Web-based single-cell databases and analytic tools, we envision that scTIME will deepen our understanding of the molecular and cellular mechanisms underlying the different responses of patients with cancer to immune therapies.

D. Hu reports grants from National Natural Science Foundation of China and Tianjin Natural Science Foundation during the conduct of the study. X. Gao reports grants from National Natural Science Foundation of China and CAMS Initiative for Innovative Medicine during the conduct of the study. No disclosures were reported by the other authors.

F. Hong: Data curation, formal analysis, validation, visualization. Q. Meng: Data curation, formal analysis, validation, visualization. W. Zhang: Data curation, formal analysis, visualization. R. Zheng: Data curation, validation. X. Li: Data curation. T. Cheng: Resources, writing–review and editing. D. Hu: Conceptualization, supervision, funding acquisition, writing–review and editing. X. Gao: Conceptualization, formal analysis, supervision, funding acquisition, investigation, visualization, writing–original draft, project administration, writing–review and editing.

The authors thank My Health Gene Co., Ltd. for their help in collecting published data and developing the Web interface for the scTIME Portal.

This work was supported by the National Natural Science Foundation of China (32000469 and 31872825), Tianjin Natural Science Foundation (18JCYBJC42400), and the CAMS Initiative for Innovative Medicine (2018-I2M-AI-010).

The publication costs of this article were defrayed in part by the payment of publication fees. Therefore, and solely to indicate this fact, this article is hereby marked “advertisement” in accordance with 18 USC section 1734.