{'fraction_non_alphanumeric': 0.026702517268555003, 'fraction_numerical': 0.011265124472671643, 'mean_word_length': 4.913377192982456, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Estimation of immune cell content in tumor using single-cell RNA-seq reference data

Background The rapid development of single-cell RNA sequencing (scRNA-seq) provides unprecedented opportunities to study the tumor ecosystem that involves a heterogeneous mixture of cell types. However, the majority of previous and current studies related to translational and molecular oncology have only focused on the bulk tumor and there is a wealth of gene expression data accumulated with matched clinical outcomes. Results In this paper, we introduce a scheme for characterizing cell compositions from bulk tumor gene expression by integrating signatures learned from scRNA-seq data. We derived the reference expression matrix to each cell type based on cell subpopulations identified in head and neck cancer dataset. Our results suggest that scRNA-Req-derived reference matrix outperforms the existing gene panel and reference matrix with respect to distinguishing immune cell subtypes. Conclusions Findings and resources created from this study enable future and secondary analysis of tumor RNA mixtures in head and neck cancer for a more accurate cellular deconvolution, and can facilitate the profiling of the immune infiltration in other solid tumors due to the expression homogeneity observed in immune cells.

Similar to the data analysis presented in Purma study, we choose to use the t-SNE method to visualize the cell clusters and explore the cell type compositions based on transcriptomes of all examined cells. However, as shown in the previous analysis and in the results section, t-SNE method alone is only able to identify clusters of major cell types and not able to distinguish between T cell subpopulations. Furthermore, the location of the clusters in the t-SNE map and their relative positions to other clusters will change across analysis runs. As a limitation of the technique, t-SNE cannot reproduce the same clustering map if different cells or perplexity parameters are chosen in one analysis run. Therefore, we propose to use a multi-stage cell identification scheme for obtaining more accurate cell type inference--by adaptively integrating t-SNE and ssGSEA results. The steps and detailed parameters used are described below.
(1) Tumor cell classification: To classify HNSCC malignant cells, we performed t-SNE analysis of all cells using perplexity parameter of 50 followed by DBscan clustering (with parameters eps=5 and minPts =5). Clusters were classified as malignant cells and nonmalignant cells based on their ssGSEA enrichment scores using signature genes for HNSCC tumor cells ( Figure S1A,B). As reported previously in various cancer studies [1,3], malignant cells were clustered by patients while non-malignant cells were clustered by cell types ( Figure S1C). (3) T cell subtype identification: Similar procedure was used to classify T cell subtypes from the lumped T cells population identified in step 2. We performed single-cell consensus clustering (SC3) analysis [8] and were able to identify four distinct clusters of T cell subpopulations. These four clusters were assigned to conventional CD4 + T cells (CD4 + Tconv), T-regulatory cells (Treg), conventional CD8 + T cells (CD8 + Tconv), and exhausted CD8 + T cells, based on their ssGSEA enrichment scores ( Figure S3A and Figure S3B). Next, differential expression analysis was performed comparing CD4 Tconv vs. Treg cells, and CD8 + Tconv vs. exhausted CD8 + T cells using R package limma [9]. Only genes with |log2FoldChange| >1 and Benjamini-Hochberg adjusted p-value < 0.05 were considered significantly differentially expressed and reported in Table S3. The identified differentially expressed genes were compared with previously reported marker genes for these cell types.

scRNA-derived marker genes
To develop a finer panel of cell-type-specific genes, we identified marker genes that are specifically expressed in each cell type. Differential expression analysis was first performed between any pairs of the 11 cell types using R package limma. Then marker genes of each cell type were identified as those significantly highly expressed in cell type under consideration compared to at least 5 other cell types (log2FoldChange >3 and Benjamini-Hochberg adjusted p-value < 0.05). In total, we identified 581 marker genes and reported the gene names and limma results in Table S4.

Deconvolution method for bulk tumor
The objective of the deconvolution algorithm is designed to solve for the linear equations derived from microarray data. In this study, we will apply the same SVR method implemented in CIBERSORT to infer cell types that are more representative in head and neck tumors. The reference GEP panels used in SVR will be described in the following section.

In silico assessment of final reference GEP panels
With the availability of high-resolution scRNA-seq data, one main objective of this study is to explore new ways to generate the reference GEP matrices to be used in bulk tumor deconvolution, i.e., the matrix B as described in the previous section. The ideal B matrix should be able to yield maximal and robust discriminatory power between cell type clusters.
Meanwhile, the pooled scRNA-seq data can be served as ground truth for benchmarking the performance of reference GEP as well as deconvolution methods-because the true cell composition in the bulk gene expression data will be known. The similar idea has been implemented in a recent study [7]. The first step of constructing reference GEP matrices is to choose a panel of reference genes that can distinguish the cell populations. In this study, we will focus on four gene panels: (1) (Table   S2). (3) The scRNA-derived marker gene panel discovered through the steps described previously in the method: which contains genes that uniquely expressed in each cell population identified from HNSC scRNA-seq data (Table S4); (4) A T-cell-specific GEP panel discovered through steps similar to GEP panel (3) but with a focus on four T cell subtypes ( Table S3). Note that we only used the gene list information of these panels. The GEP matrix of these genes is formed through averaging all single cells assigned to these populations. In order to assess the prediction performance of the above four GEP panels, we tested them on in silico bulk tumors by aggregating the single cell transcriptome data.
Expression data of individual cells from the same patient in Puram study were pooled to form 15 in-silico tumors, which exhibit varied cellular compositions. . However, results from scRNA-seq data suggests that the overall Treg expression signature may be underrepresented in genomic projects that are biased towards tumors with higher purity, such as TCGA. In the following, we briefly describe results generated from each step. First, we observed that the unsupervised clustering on all cells based on t-SNE revealed eight major clusters as depicted in Figure S1A. Note that, at this stage, we had no information about cell types underlying these cell groups and the number of clusters might differ subject to the perplexity parameter choice in t-SNE. We started the cell type identification from first distinguishing tumor and non-tumor cells. By adding ssGSEA scores representing the tumor cell signature into the t-SNE map ( Figure S1B), we identified two major cluster regions of malignant cells located in the very top and lower regions. By further adding the color layers reflecting the tumor origin, we observed that the cell clusters in these regions were clearly separated by patient IDs while they were mixed together in a mosaic pattern in other cluster regions ( Figure S1C). The results above align with previous findings [1,3,7] that inter-tumor heterogeneity may arise more at the tumor malignant cell level than at the immune cell level-suggesting that immune cell signatures abstracted from the proposed scheme will be applicable to not only HNSCC samples generated from different studies but also samples from different tumor types. Next, we performed a second round of t-SNE analysis by excluding all tumor cells identified from previous steps. The new clustering analysis revealed seven major cell clusters ( Figure S2A). We were able to identify the cell types corresponding to each cluster by adding ssGSEA score specific to Fibroblasts, B cell, Macrophages, endothelial cells, dendritic cells, mast cell, and T cells one at each time as depicted in Figure S2B. As expected, this subset of cell population is dominated by

Identifiable
Fibroblasts and T cells. When we adding the color layers reflecting patient origins into Figure   S2A, we found a similar pattern that patient IDs were mixed together in each cell type cluster, indicating that the sub-clusters (such as in the T cells) may reveal further cell subtypes. This leads us to the next step by further zooming into the expression profiles of cells from T cell populations.

Deconvolution of T cell subtypes using identified T cell population. Based on SC3, we further identified four clusters from T cells (Figures 2A). The cell types in the T cell
subpopulation were first determined based on the gene enrichment signatures of CD4 + and CD8 + cells (as shown in the upper panel in Figure S3B). Within these two subpopulations, [CD8 + cells further marked with ssGSEA signatures for CD8 + Tconv and CD8 + exhausted; and CD4 + cells were marked with CD4 + Tconv and Treg cells signature values ( Figure S3B).
As shown in Figure S3, the signatures for two CD8 + cell types are overlapped and it is difficult to assign these cells to any subtypes. As further summarized in the heatmap of ssGSEA scores ( Figure S4), the ssGSEA analysis based on curated signature genes were able to distinguish between major cell types using single cell level expression data but failed to provide the necessary granularity in separating T cell subtypes. To determine T cell subtypes, especially CD8 + subtypes, we performed differential expression analysis between the two cell groups identified within CD4 + T cells and CD8 + T cells. Differentially expressed genes (adjusted p value < 0.05, limma moderated t-test, and |log2fold-change| >1) are reported in Table S3. Cell subtypes were then inferred from the status of top differentially expressed genes, by comparing them with existing cell-type-specific marker genes. Figure   2B and Figure 2C are heatmaps depicting top differentially expressed genes between CD8 + cell clusters and CD4 + cell clusters, respectively. Candidate genes that overlapped with marker genes identified from previous studies are listed and labeled in heatmaps. Note that several exhaustion-related genes can serve as markers for separating both subtypes in CD4 and CD8, such as TIGIT and CTLA4. For the CD8 + T cell subtypes, we compared the candidate marker genes identified in our DE analysis to the exhausted CD8 + T cells marker genes reported in a previous single-cell RNA-seq from infiltrating T cells of lung cancer [13].
A total of 36 genes are found shared by the two studies and all labeled in Figure 2B. Among these 36 genes also includes 14 known exhaustion markers, such as PCCD1, TIGIT, HAVCR2, and CTLA4 ( Figure 2B, text in red), which further confirmed the identify of these exhausted CD8 + T cells. The other CD8 + T cell cluster without expression of exhaustion genes is considered as conventional CD8 + T cells. For the CD4 + T cell subtypes, we also compared the candidate marker genes identified from the DE analysis with the Tregs marker genes reported by four previously published scRNA-seq data from different cancer types [13-16] ( Figure 2D). We observed that there were 20 genes shared by all five studies ( Figure   2C, text in red), including known Tregs markers FOXP3, TIGIT, and CLTA4; and there were many more genes previously identified at least once ( Figure 2E). Our study also identified 207 genes that uniquely enriched in this HNSCC dataset ( Figure 2E), including PPP1CA, RUNX3, CCR6, and PSMB8 which were previously reported to be associated with Tregs and their functions [17][18][19][20]. Based on these observations, we assigned Tregs to this cluster of CD4 + T cells. The other CD4 + cluster with low expression of exhaustion markers and with exclusively high expression of CCR7, CXCR4, and TOBI was considered as conventional CD4 + T cells.
Evaluation of prediction performance of reference GEPs. For each cell type identified from previous steps, we established cell-type-specific reference GEP matrix by the mean expression values of selected genes. We use C1 to denote the curated gene list from previous literatures which are used in ssGSEA (Table S2), C2 to denote marker genes selected from the DE analysis described above (Table S4), T1 to denote the marker genes selected from DE analyses for separating T cell subtypes (Table S3), and M1 to denote marker genes selected from DE analyses for separating tumor and non-tumor cells. In our analysis, we constructed reference GEP matrices by taking the mean from the following ensemble gene lists: (1) LM22, (2) C1, (3) C2, (4) LM22+C1, (5) LM22+C1+T1, (6) LM22+C1+T1+M1, and (7) LM22+C1+C2+T1+M1. As presented in Figure S5, we evaluated the prediction performance of CIBERSORT using these GEPs in terms of correlation between predicted abundance and the true abundance in the simulated bulk tumor (through pooling all cells in one patient, see Methods). We observed that all of these reference GEPs achieved promising prediction accuracies (r>0.9). This result indicates that existing marker genes provides saturated signatures if forming GEPs on right cell groups. Therefore, we will focus on the evaluation of the LM22+C1 gene panel because of it has a moderate number of genes and all genes included are well studied. All reference GEPs matrices used in this study are provided in Supplementary Table S5.
Scatterplots in Figure 3A demonstrate strong correlations between true cell proportions and predicted cell proportions based on GEP curated form LM22+C1 scRNA-seq data, where each point represents a simulated bulk sample. Figure 3B further compares the cell abundance estimation accuracy (correlation) for the reference GEP included in CIBERSORT and the reference GEP trained based on the LM22+C1 scRNA-seq panel. Our method shows better prediction performance in all case for cell types that CIBERSORT can provide estimation, especially in estimating CD8 T cells. We further gauged the estimated cell proportion from CIBERSROT by taking into account the fact that the original GEP only include reference for immune cells. Such adjustment was made by assuming that tumor cell (purity) and stromal cell proportion were known so that a relative abundance on each remaining cell types can be calculated. Even with this unrealistic scenario, the prediction performance based on the adjusted proportion was still inferior to the scRNA-seq trained GEP in all cases. But we did observe that CIBERSORT estimation on macrophages and dendritic cells was greatly improved with this adjustment ( Figure S6). To test the robustness of the GEP panel to the cell components, we re-run the devolution analysis on all simulated samples using the leave-one-out GEP, i.e. each time we remove one cell-type-specific vector from the GEP matrix. As shown in Figure S7, the high prediction accuracy was maintained in most scenarios, and only the estimations for fibroblasts and malignant cells were detectably impacted by the leave-one-out GEP.
Although C2 and T1 gene sets (determined based on DE tests) did not provide additional information as a gene panel in constructing GEP, they provide a new alternative cell-typespecific biomarker for future studies. As shown in violin plots ( Figure S14), these markers are exclusively over-expressed in cell types that they are representing, indicating their validity as independent surrogate biomarkers. A total of 182 genes were found overlapping between groups C2+T1 and LM22+C1. Expressions of these genes for each single cell were plotted in Figure S15, demonstrating their ability as biomarker panel alone to separate major cell types but not T cell subtypes.
Finally, as a supportive validation, we tested the proposed scGEP on TCGA HNSCC tumor samples and compared with results generated from similar methods developed for bulk tumor deconvolution. Figure S8 compared  proportion estimates from the scGEP-based method. As shown in Figure S9, the analysis showed a good agreement between two methods. We also compared the estimated total immune cell proportions and total T cell proportions between HPV positive and HPV negative cancer patients. As expected, tumors from HPV positive patients showed higher infiltration of immune cells and T cells ( Figure S10). Abundance of tumor infiltrating CD8 and total immune cells were also found associated with survival outcomes in TCGA HNSCC patients ( Figure S11).  Figure S12), but had limited impact on other cell types ( Figure S13).

Discussion
One notable observation in our simulation studies is that the GEP calculated based on existing marker gene panel (LM22+C1) can provide as accurate a predictive capability as the genes selected only from the differential expression of single cell populations, although they are overlapped in many genes. We conclude that the prediction performance of GEP is more sensitive to the cell populations purified for a particular cell type than the marker gene panel. Nevertheless, the newly discovered genes from the scRNA-seq data and their underlying pathway warrant further validations as potential biomarkers, especially those genes that are differentially expressed between T subtypes.
Although we have only tested support vector regression method for cell mixture estimation, the HNSCC single cell sequencing data curated from this study provides a useful source for the assessment of accuracy of newly developed deconvolution methods. For example, the core SVR algorithm implemented in CIBESRORT only uses a single kernel under the fixed default parameter setting. The prediction performance might be improved through searching for an optimal kernel or using the state-of-the-art multiple kernel learning technologies [24]. Currently there is a lack of suitable benchmark dataset that allows a fair and systematic evaluation of methods for estimating cell mixtures in solid tumors. Weak correlations were often found between molecular-data-based estimations and pathology based methods such as IHC and H&E images [29]. This is partially due to the fact that each of these assays was carried out using input materials from different parts of a tumor.
Because all cell proportions are known, the in silico pooled bulk tumor data from individual cells provides a more accurate reference at almost zero cost. Plus the composite cells from a single tumor could better mimic the real case scenario than creating bulk expression dataset through conducting RNA-seq on randomly mixed cells. For head and neck cancer per se, the scRNA-seq data from Puram study provide an ideal source for both training and validation purposes because the studied tumors have (1) uniformly varied tumor purity, and (2) it provides reference for subpopulations such as exhaustive CD8 + T cell that were not present in previous scRNA-seq experiments on melanoma and lung cancers. A limitation of the in silico method is that the cell size factor has not been taken into account. As cell types of different size have different amount of RNA yield, it is of interest for future research to be able to adjust for the cell size factor so that the estimated relative abundance will be closer to absolute cell proportions.
The key idea proposed in this work is most similar to a previous study conducted by Schelker et al. [7] which focused on scRNA-seq data from melanoma and PBMC. The two main differences between the two works are (1) the melanoma data used by Schelker et al.
only provided sufficient information for distinguish nine major cell types and three T cell subtypes, whereas the HNSCC data we studied was able to further separate exhaustive CD8 T cells and provide corresponding reference GEP; (2) In our method, we used both marker gene information and a global ssGSEA scores to determine cell types from adaptive clustering analysis. We believe that more studies along this line will be conducted to generate more accurate cancer-type-specific and T-cell-subtype-specific reference GEP.
Finally, we believe that apart from looking for reference profiles based on gene expression, the same approach can be extended in future search to identify reference DNA methylation profiles (DMP). DMP will be a promising new resource for tumor composition deconvolution because Alternations at DNA methylation level are deemed to be more stable than the gene expression level. But the single-cell DNA methylation analysis, such as bisulfite sequencing, is still in an experimental phase.

Ethics approval and consent to participate
Not applicable.

Consent for publication
Not applicable.

Availability of data and materials
The HNSCC scRNA_seq data used in this study was retrieved from Puram et al [1].

Competing interests
The authors declare that they have no competing interests.

Funding
This work was supported in part by Institutional Research Grant number 14-189-19 from the American Cancer Society, and a Department Pilot Project Award from Moffitt Cancer Center.
The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.