Single-cell RNA sequencing demonstrates the intratumoral heterogeneityof papillary thyroid carcinoma

Sample collection and clinical information

Four patients diagnosed with Papillary thyroid carcinoma were obtained from Shanghai Tenth People’s Hospital(Supplementary Data 1). The average age was 37.5 years old and 50% of them were female. A total of 15 tissue samples (7 tumor tissue, 6 distant normal tissue, and 2 metastatic lymph node) were obtained from 4 Papillary thyroid carcinoma patients.

Ethics statement

This study has complied with all relevant ethical regulations for the human participants.

Tissue digestion and single cell suspension preparation

Cells of each sample were firstly stained with two fluorescent dye, Calcein AM (Thermo Fisher Scientific Cat. No. C1430) and Draq7 (Cat. No. 564904), for precisely determination of cell concentration and viability via BD Rhapsody™ Scanner before single cell multiplexing labelling. Cell of each sample were sequentially labeled with BD Human Single-Cell Multiplexing Kit (Cat. No. 633781) which utilizing an innovative antibody-oligo technology ¹⁴ mainly to provide higher sample throughput and eliminate batch effect for single cell library preparation and sequencing. A set of 12 antibodies in this kit recognize the same universally expressed cell-surface antigen of human cells. Each antibody is conjugated with a Sample Tag, a unique 45-nucleotide barcode sequence. Briefly, cells from each sample were labeled by antibodies with different sample tags respectively and washed twice with BD Pharmingen™ Stain Buffer (FBS) (Cat. No. 554656) before pooling all samples together. BD Rhapsody Express system¹⁵based on Fan et al was utilized for single cell transcriptome capturation. Pooled samples were then loaded in one BD Rhapsody™ Cartridge that was primed and treated strictly following the manufacturers protocol (BD Biosciences). Cell Capture Beads (BD Biosciences) were then loaded excessively onto the cartridge to ensure that nearly every micro-well contains one bead, and excess beads were washed away from the cartridge. Viable cells with beads captured in wells were detected in BD Rhapsody™ Scanner and the total number was 10871 and cell multiplet rate was 3.00%. Then cells were lysed, Cell Capture Beads were retrieved and washed prior to performing reverse transcription and treatment with exonuclease I.

Library construction and sequencing

Transcriptome and SampleTag information of single cells was obtained through BD Rhapsody System. Microbeads-captured single cell transcriptome and SampleTag sequences were generated into cDNA library and SampleTag library separately containing cell labels and unique molecular identifiers (UMI) information. Briefly, double strand cDNA was firstly generated from microbeads-captured single cell transcriptome through several steps including reverse transcription, second strand synthesis, end preparation, adapter ligation and whole transcriptome amplification. Then final cDNA library was generated from double strand full length cDNA by random priming amplification with BD Rhapsody cDNA Kit (BD Biosciences, Catalog No.: 633773) and BD Rhapsody Targeted mRNA & AbSeq Amplification Kit (BD Biosciences, Catalog No.: 633774). On the other hand, SampleTag library was generated from microbeads-captured single cell SampleTag sequences through several steps including reverse transcription, nest PCR and final index PCR. All the libraries were sequenced in a PE150 mode (Pair-End for 150bp read) on the NovaSeq platform (Illumina).

Raw data preprocessing and quality control

Raw sequencing reads of cDNA library and SampleTag library were processed through the BD Rhapsody Whole Transcriptome Assay Analysis Pipeline (Early access), which included filtering by reads quality, annotating reads, annotating molecules, determining putative cells and generating single-cell expression matrix. Briefly, read pairs with low sequencing quality were firstly removed. The quality-filtered R1 reads were analyzed to identify sequence of cell labels, UMI sequence, and poly-dT tail sequence, meanwhile the quality-filtered R2 reads were mapped to Genome Reference Consortium Human Build 38 (GRCh38) using STAR (version 2.5.2b) in reads annotation step. Further adjustment was performed with the recursive substitution error correction (RSEC) and distribution-based error correction (DBEC) algorithms to remove artifact molecules due to amplification bias in molecules annotation step. Putative cells were distinguished from background noise with second derivative analysis in putative cell determination step. Finally, putative cells information was combined together with molecules adjusted by the recursive substitution error correction and distribution-based error correction algorithms to generate a single-cell expression matrix. The pipeline also determined sample origin of every single cells via Sample determination algorithm according to the sequencing reads of SampleTag library. The algorithm classified all putative single cell into three categories: Multiplet, two or more SampleTag exceed their minimum thresholds, indicating more than one actual cell in one micro-well. Undetermined, not enough SampleTag reads for the sample origin. SampleTag01-12, only one SampleTag exceed their minimum thresholds.. Among all the output files, Matrix for UMI counts per cell corrected by DBEC and SampleTag annotation result were used for downstream clustering analysis.

Dimension reduction and clustering

The R package Seurat¹⁵ was utilized for subsequent analysis. Raw gene expression matrices from the cartridge were read into R (version 3.6.0) and converted to Seurat objects. Cells label as “Undetermined” and “Multiple” were excluded in the following analysis. The gene expression matrix was then normalized to the total cellular UMI count. Top 2000 features were selected as highly variable genes for further clustering analysis. In order to reduce dimensionality, PCA was performed based on the highly variable genes after scaling the data with respect to UMI counts. On top of that, the first 50 principle components were chosen for downstream clustering based on heatmap of principle components, Jackstraw plot, and elbow plot of principle components to further reduce dimensionality using the tSNE algorithm. The transcriptional markers of each cluster were calculated using the FindAllMarkers function with the Wilcoxon Rank-Sum Test under the following criteria: log2 fold change > 0.25; p < 0.01; min. percentage > 0.25. Top markers of each cluster were then selected to perform heatmap plot.

Cell annotation and cell type identification

Cell populations were matched to cell types based on the expression of known marker genes and previously identified expression signatures^{15, 16, 17, 18, 19, 20}.

Comparison of cell clusters and cell type proportion

The change in fraction of the different cell populations (clusters) was separately computed for each sample across all clusters, as the fraction of cell in each cluster, out of the total number of cell²¹. To assess statistically significant changes in a fraction of a specific population, we first performed the Shapiro–Wilk test for normality (using the shapiro.test function in R), followed by paired t-test (for normally distributed populations) and an additional paired nonparametric Wilcoxon test (using t.test and wilcox.test functions in R).

Significantly dysregulated genes identification

We identified differentially expressed genes (DEGs) based on analysis of variance (ANOVA) using the R. In brief, only genes that met these criteria were considered as DEGs: 1)P-value of F test < 0.05; 2) |FC| ≥2.

External data validation with TCGA datasets

CIBERSORTx is a new machine-learning method designed for detecting the abundance of cell types in bulk RNA-seq data; this method has been called “digital cytometry” by the authors²². To evaluate the roles of the cell clusters identified in the present work, we used CIBERSORTx to analyze TCGA expression matrix, which was first normalized to obtain log2(TPM+1) values. As the signature matrix was constructed, we calculated the relative abundance of each cell cluster and then divided the patients into two groups: high 50% and low 50%. Subsequently, Kaplan-Meier analysis was performed using ggsurvplot function of the R package ‘survminer’. One-way ANOVA was used to evaluate whether the occurrence of cell clusters in the different subgroups showed significant differences.

Functional annotation and enrichment

GO enrichment and KEGG enrichment of cluster markers were performed using KOBAS software with Benjamini-Hochberg multiple testing adjustment, using the top 20 markers gene of cluster. The results were visualized using R package.

Gene set variation analysis (GSVA) ²³was performed using 50 hallmark gene sets obtained from the molecular signature database using default sets, as described in the GSVA package.

Identification of hub regulons with SENIC

In order to investigate the gene regulatory network of the samples, we utilized by SCENIC²⁴, a python-based algorithm for the reconstruction of the transcriptional states and regulatory networks from the scRNA-seq data,37 in the Jupyter notebook software (version1.0.0).

Cell-cell communication analysis based on ligand-receptor pairs

CellPhoneDB²⁵ is a Python-based computational analysis tool; it enables the analysis of cell-cell communication at the molecular level. A website version was also provided for the analysis of relatively small datasets (https://www.cellphonedb.org/).

Cell trajectory analysis with Monocle

In order to reveal the changes of immune cells in the tumor-educating process, we employed Monocle 2, an R package designed for singlecell trajectories²⁶. The following parameters were set: mean expression R0.125, num_cells_expressed R10, qval <0.01 (differentialGeneTest function). The trajectories were visualized as 2D tSNE plots, while dynamical expression heatmaps were constructed using the plot_pseudotime_heatmap function.

Drug-target analysis with TCMID database

After obtaining the differentially expressed genes, we obtain the targeted relationship between drugs and genes from the TCMID database(http://www.megabionet.org/tcmid).

CNV Estimation

CopyKAT²⁷(Copynumber Karyotyping of Tumors) is a computational tool using integrative Bayesian approaches to identify genome-wide aneuploidy at 5MB resolution in single cells to separate tumor cells from normal cells, and tumor subclones using high-throughput sc-RNAseq data. Cells with extensive genome-wide copy number aberrations (aneuploidy) are considered as tumor cells, wherease stromal normal cells and immune cells often have 2N diploid or near-diploid copy number profiles. In this study, we inferred tumor subclones from using the R package.

Statistical analysis

Box plots were generated using the R base package and default parameters. Hence, the boxes span the interquartile range (IQR; from the 25th to the 75thpercentiles), with the centerline corresponding to the median. Lower whiskers represent the data minimum or the 25th percentile minus 1.5×IQR, whichever is greater. Upper whiskers represent the data maximum or the 75th percentile plus 1.5×IQR (lower), whichever is lower. Violin plots were generated using the beanplot R package, and data distribution band width was estimated by kernel density estimation, as per the built-in ‘nrd0’ option. Comparisons between two groups were done using unpaired two-tailed t-tests. One-way analysis of variance (ANOVA) with Tukey’s multiple comparisons tests were used for multiple group comparisons.

Single cell atlas and heterogeneity of papillary thyroid carcinoma (PTC)

A total of 15 samples from four papillary thyroid carcinoma patients were involved in this study. After quality filtering, 30986 genes were detected, of these, 13390 cells originated from 7 primary papillary thyroid carcinoma, 2869 cells from 2 lymph metastatic, and 11945 cells from 6 paired normal thyroid tissue samples (Fig. 1A,B). Subsequently, we partitioned the cells into 26 clusters, and classified those clusters into 10 major cell types based on known markers described in previous studies: B cells (CD19+MS4A1+CD38+CD79A+CD79B+); CD4+T cells (CD3D+CD4+); CD8+T cells

• Download figure
• Open in new tab

Figure 1. Overall design and single cell atlas in thyroid carcinoma.

A Workflow diagram showing the processing of samples. B UMAP plots of total cells, colored by sample origin (primary papillary thyroid carcinoma, lymph metastatic or paired normal thyroid), the corresponding patient, the associated cell type. C Expression of cell marker genes. Additional cell marker genes were shown in Supplementary Fig. 1. D Changes in frequency of multiple cell types and clusters in primary papillary thyroid carcinoma and lymph metastatic. For each cell type and cluster, the average fraction of cells from each patient group is shown, Asterisks on the left of the vertical line denote statistically significant differences between primary papillary thyroid carcinoma and normal thyroid tissue, meanwhile, asterisks on the right is used to show statistically significant differences between lymph metastatic andprimary papillary thyroid carcinoma. *p < 0.05, **p < 0.01, ***p < 0.001, two-tailed t-tests. E Interaction network constructed by CellPhoneDB; size of circles and color of arrows represent interaction counts, brighter color and larger size means more interaction with other cell types.

(CD3D+CD8A+); endothelial cells (CD31+CD34+); epithelial cells (EPCAM+KRT18+); fibroblasts (COL1A1+); myeloid cells (CD14+CD86+ITGAX+CD80+CD83+ITGAM+); naive T cells (CD3D+CCR7+); NKT cells (CD3D+NKG7+); and plasma cells (CD79A+SDC1+) (Fig. 1C, Supplementary Fig. 1).

By calculating the frequency of the cell types, we revealed similar cellular landscapes in primary papillary thyroid carcinoma and normal thyroid tissue samples, with differences in the proportions of several cell types (P<0.01; Fig. 1D). There were major differences in immune cells in primary tumor relative to normal tissue, the frequency of NKT cells/plasma cells increased in primary tumor (P<0.01; Fig. 1D), while that of CD4+ Tcells/B cells decreased. The most striking changes between lymph metastatic and primary tumor were observed in CD8+T cells. At the cluster level, there were also differences in the frequency of clusters within the same cell group, take myeloid cells as an example, we observed an increase in the frequency of cluster 21 during the disease progression, at the same time, decrease in the proportion of cluster 25 was also be identified. These results demonstrate that the disease progression is associated with major changes in immune cells,heterogeneity is common between samples and cell populations.

To investigate the interaction network in the tumor microenvironment in papillary thyroid carcinoma, CellphoneDB was be used to calculate potential ligand-receptor pairs in Cells. Cytoscape was performed to visualized the cell interaction. We found that myeloid cells and T cell related cells possessed the most interaction pairs with cells from other lineages (Fig. 1E), demonstrated the dominant role of myeloid cells and T cells.

Identification of malignant cells in epithelial cells

In order to distinguished malignant from non-malignant cells within the epithelial cells, CopyKAT was performed to identified papillary thyroid carcinoma genome alterations. CNVs accumulated in most lymph metastatic papillary thyroid carcinoma,despite the heterogeneity, almost all malignant cells possessed deletions from chromosomes 16 and 19 and amplifications in chromosomes 13(Fig. 2A,B). We distinguished 1679 nonmalignant cells, and 450 cells were identified as malignant cells, of which malignant cells were mainly composed of lymphatic metastasis samples, but we also found very few malignant cells in paired normal thyroid tissue samples(Fig. 2B). In order to support the identification results of the malignant cells, pseudotime trajectory analysis was performed. The results showed that malignant cells was present at the end of the differentiation trajectory(Fig. 2C).

• Download figure
• Open in new tab

Figure 2. Identification of malignant cells in epithelial cells.

A Clonal substructure of epithelial cells delineated by clustering single-cell copy number profiles inferred from scRNA-seq data by CopyKAT. B UMAP plots of total epithelial cells, colored by the associated cell type, sample origin (primary papillary thyroid carcinoma, lymph metastatic or paired normal thyroid), the frequency of cell types and clusters in malignant and non-malignant cells. C Differentiation trajectory ofepithelial cells, with each color coded for pseudotime (right) and clusters (left). D GO enrichment analysis of up-regulated genes in malignant lymphatic metastasis cells. while color represents expression level. E Heatmap of the area under the curve (AUC) scores of TF motifs estimated per cell by SCENIC. Shown are differentially activated motifs in malignant and non-malignant cells, respectively.

We further characterized the functions of differential genes between lymphatic metastasis and non-metastasis in malignant cells by comparing pathway activities, pathways involved in cellular component assembly, protein containing complex assembly and regulation of microtubule motor activity were relatively upregulated in lymphatic metastasis malignant cells(Fig. 2D).In order to screen for key genes related malignant cells,we conducted SCENIC analysis,Via SCENIC analysis, we identified ATF3 and BHLHE40 motifs were highly activated in malignant cells(Fig. 2E).These results provide potential targets for suppressing cells to develop into malignant cells.

Increased CD4+Tregs in tumor and lymph metastatic PTC and decreased transitional EGR1+CD4+T cells

For discussing the heterogeneity among T cells, 17690 T cells were clustered into 10 subclusters based on known markers described in previous studies: CD4+Tregs(CD3D+CD3E+CD3G+CD4+FOXP3+); Cytotoxic CD8+T cells(CD3D+CD3E+CD3G+CD8A+NKG7+); Exhausted CD4+T cells(CD3D+CD3E+CD3G+CD4+PDCD1+CTLA4+TIGIT+LAG3+); Follicular helper (Tfh) T cels(CD3D+CD3E+CD3G+ICOS+); Naïve CD4+T cells(CD3D+CD3E+CD3G+CD4+CCR7+); Naïve T cells(CD3D+CD3E+CD3G+CCR7+);

NKT cells(CD3D+CD3E+CD3G+NKG7+GZMB+GNLY+) ; Pre-exhausted CD8+ T cells(CD3D+CD3E+CD3G+CD8A+PDCD1+CTLA4+TIGIT+LAG3+); Proliferating CD8+T cells(CD3D+CD3E+CD3G+CD8A+MKI67+TOP2A+CDK1+PCNA+); Transitional EGR1+CD4+T cells(CD3D+CD3E+CD3G+CD8A+MKI67+TOP2A+CDK1+PCNA+)(Fig. 3A,B). In addition to the marker used above, we also found some other makers that can be used for cell population identification In the future, for example, TYMS can be used for Pre-exhausted CD8+ T cells identification, RTKN2 can be used for Pre-exhausted CD4+Tregs identification(Fig. 3C). Through cell frequency analysis, we found an interesting phenomenon, the proportion of CD4+Tregs continues to increase during the development of the disease, while the trend of Transitional EGR1+CD4+T cells is the opposite(Fig. 3D). The trend was consistent with TCGA data calculated by cibersortx(Fig. 3E). The results indicate that both of CD4+Tregs and Transitional EGR1+CD4+T cells may have an important impact on disease progression, and our newly identified cell population(Transitional EGR1+CD4+T cells) may be antagonistic to the CD4+Tregs.

• Download figure
• Open in new tab

Figure 3. T cell population heterogeneity.

A UMAP plot of T cells color-coded by their associated clusters. B Dotplot of cell markers; sizes of dots represent abundance, while color represents expression level. C Potential new cell markers, additional new cell marker genes were shown in Supplementary Fig. 2. D Changes in frequency of multiple cell types and clusters in T cell population. For each cell type and cluster, the average fraction of cells from each patient group is shown, Asterisks on the left of the vertical line denote statistically significant differences between primary papillary thyroid carcinoma and normal thyroid tissue, meanwhile, asterisks on the right is used to show statistically significant differences between lymph metastatic andprimary papillary thyroid carcinoma. *p < 0.05, **p < 0.01, ***p < 0.001, two-tailed t-tests. E Violin of cell abundance predicted from TCGA THCA cohort by CIBERSORTx.

As shown in Fig. 4A, the differentially expressed genes derived from transitional EGR1+CD4+T cells could be identified, 31 mRNAs were found to be up-regulated and 50 were down-regulated in primary tumor compared to normal thyroid tissue samples(Fig. 4A). We further characterized the functions of CD4+Tregs and Transitional EGR1+CD4+T cells by comparing pathway activities. in primary tumor, CD4+Tregs and Transitional EGR1+CD4+T cells exhibited large different activated pathways, pathways involved in epithelial-mesenchymal transition and hypoxia were relatively upregulated in the CD4+Tregs, whereas G2M checkpoint was down regulated in Transitional EGR1+CD4+T cells(Fig. 4B). In lymph metastatic, mitotic spindle and oxidative phosphorylation were also down regulated in Transitional EGR1+CD4+T cells(Supplementary Fig. 4).In order to screen for key genes related to tumor occurrence and development in CD4+Tregs and transitional EGR1+CD4+T cells,we conducted SCENIC analysis,Via SCENIC analysis, we identified essential motifs in CD4+Tregs and transitional EGR1+CD4+T cells. CREM, FOSL2, ATF3, REL, and HES1 are CD4+Tregs-specific motifs. BCLAF1 and ETS1 motifs were highly activated in Transitional EGR1+CD4+T cells(Fig. 4C).These results provide potential targets for inhibiting or reversing the formation of the immunosup-pressive microenvironment.

• Download figure
• Open in new tab

Figure 4. CD4+Tregs and transitional EGR1+CD4+T cells in tumor and lymph metastatic PTC and potential therapeutic target.

A Volcano plot of differentially expressed genes (DEGs) between lymph metastatic and primary tumor in transitional EGR1+CD4+T cells. Upregulated genes and downregulated genes were colored in red and blue, respectively. FC≥2. Additional Volcano plot were shown in Supplementary Fig. 3. B Differences in 50 hallmark pathway activities scored with GSVA software. Shown are t values calculated by a linear model. C Heatmap of the area under the curve (AUC) scores of TF motifs estimated per cell by SCENIC. Shown are differentially activated motifs in CD4+Tregs and transitional EGR1+CD4+T cells, respectively. D Targeted drug prediction corresponding to differential genes in TCMID database, Additional Drug-Target Genes were shown in Supplementary Fig 5. E Heatmap of positive and negative immune checkpoint expression on T cells in lymph metastatic PTC. Additional immune checkpoint expression were shown in Supplementary Fig 5

We extracted small molecule drugs that can target differential genes related to Transitional EGR1+CD4+T cells from the TCMID database(Fig. 4D). Subsequently, we examined the immune checkpoints in the cell clusters (Fig. 4E). Notably, in lymph metastatic, we found that inhibitory checkpoints,Notably, we found that inhibitory checkpoints, TIGIT was upregulated in cells from Exhausted CD4+T cells, CD96 was upregulated in cells from Pre−exhausted CD8+T cells, since these molecules are markers of T cell exhaustion, these data indicated that cells from Exhausted CD4+T cells and Pre− exhausted CD8+T cells were exhausted in the tumor microenvironment. This result is consistent with our previous cell population classification results. Since the expression of TIGIT and CD96 was much higher than that of PDCD1 in exhausted T cell types, our data indicated that TIGIT and CD96 may be better targets for immune therapy in lymph metastatic papillary thyroid carcinoma. Opposite of this, PDCD1 maybe the most suitable targets of immunotherapy in primary papillary thyroid carcinoma.

Neutrophils cells are extremely reduced in the tumor

For discussing the heterogeneity among myeloid cells, 1607 T cells were clustered into 5 subclusters based on known markers described in previous studies:Macrophage(MSR1+); Dendritic cells(CD83+); monocyte(FCGR3A+); neutrophil(S100A8+S100A9+IFITM2+FCGR3B+CXCR1+CXCR2+); myeloid-derived suppressor cells(HLA.DRB5+) (Fig. 5A,B). Through cell frequency analysis, we found decreased proportion of neutrophils in tumors, however the proportion of macrophage continues to increase during the development of the disease(Fig. 5C), The trend was consistent with TCGA data calculated by cibersortx(Fig. 5D). We also found the DFS of the low-dendritic cells group was longer than that of the high-Dendritic cells group in the TCGA THCA(Fig. 5E), the occurrence of dendritic cells was related to bad prognosis in the case of patients from TCGA THCA cohort. The results suggested that macrophage, neutrophils, and Dendritic cells may be more relevant to papillary thyroid carcinoma.

• Download figure
• Open in new tab

Figure 5. Myeloid cells population heterogeneity.

A UMAP plot of Myeloid cells color-coded by their associated clusters. B Dotplot of cell markers; sizes of dots represent abundance, while color represents expression level. C Changes in frequency of multiple cell types and clusters in Myeloid cell population. For each cell type and cluster, the average fraction of cells from each patient group is shown, Asterisks on the left of the vertical line denote statistically significant differences between primary papillary thyroid carcinoma and normal thyroid tissue, meanwhile, asterisks on the right is used to show statistically significant differences between lymph metastatic andprimary papillary thyroid carcinoma. *p < 0.05, **p < 0.01, ***p < 0.001, two-tailed t-tests. D Violin of cell abundance predicted from TCGA THCA cohort by CIBERSORTx. E Kaplan-Meier curves of DFS for patients with THCA based on the percentage of dendritic cells in the TCGA datebase.

We further characterized the functions of myeloid cells sub types by comparing pathway activities. in primary tumor, pathways involved in mtorc1 signaling and E2F targets were relatively upregulated in the myeloid-derived suppressor cells, whereas The hallmark related pathways in neutrophils were generally down-regulated(Fig. 6A), which is consistent with the changes in cell proportions we observed earlier(Fig. 5C). However, in lymphatic metastasis samples, myeloid-derived suppressor cells were down-regulated in some hallmark related pathways(Fig. 6B). The results show that in the process of disease progression, myeloid cells are generally up-regulated in hallmark related pathways, but there is a difference between lymphatic metastasis and primary cancer.In order to screen for key genes related to tumor occurrence and development in myeloid cells,we conducted SCENIC analysis,by means of SCENIC analysis, we identified essential motifs in myeloid cells sub types. The overall regulatory factors in neutrophil were also inconsistent with other cell types(Fig. 6C). we recognized that the activity of a lot of motifs, such as YY1 and IRF5, were downregulated, while activation of the STAT3, ZNF143, and HCFC1 motifs led to the reduction of neutrophil.

• Download figure
• Open in new tab

Figure 6. M2 macrophages are strongly enriched in the tumor.

A,B Heatmap shows difference in pathway activities scored by GSVA per cell between different sample groups. Shown are t values calculated by a linear model. C Heatmap of the area under the curve (AUC) scores of TF motifs estimated per cell by SCENIC. Shown are differentially activated motifs in myeloid cells sub types, respectively. D Dotplot of M1 and M2 macrophages markers; sizes of dots represent abundance, while color represents expression level.

As shown in Fig 6D, cells from cluster 8 uniquely expressed the M1 macrophages marker CCL5, whereas M2 macrophages markers, including CTSB, CTSD, and FN1, were expressed in cluster 0.The proportion of cluster 0 is high in tumor and lymph samples, although the difference has not reached a significant level(Fig. 5C). All of this evidence certified that cluster 0 represented an M2-like TAM cluster, the increase in the proportion of cell subtypes may be related to disease progression. Cells from cluster 9 highly express the cytokine CCL22(Supplementary Fig. 6), cytokine CCL22 binds to CCR4 which expressed on the cell membrane of Tregs and has strong chemotaxis to Tregs. The results indicate that the dendritic cell cluster 9 can recruit tregs into the tumor area, thereby acting as immunosuppressive effect.