Immunologic Signatures across Molecular Subtypes and Potential Biomarkers for Sub-Stratification in Endometrial Cancer

Current molecular classification approaches for endometrial cancer (EC) often employ multiple testing platforms. Some subtypes still lack univocal prognostic significance, highlighting the need for risk sub-stratification. The tumor immune microenvironment (TIME) is associated with tumor progression and prognosis. We sought to investigate the feasibility of classifying EC via DNA sequencing and interrogate immunologic signatures and prognostic markers across and within subtypes, respectively. Formalin-fixed paraffin-embedding (FFPE) samples from 50 EC patients underwent targeted DNA and RNA sequencing, and multiplex immunofluorescence assay for TIME. DNA sequencing classified 10%, 20%, 52%, and 18% of patients into the subtype of POLE-mutant, microsatellite instability-high (MSI-H), TP53-wt, and TP53-mutant. POLE-mutant tumors expressed the highest T-effector and IFN-γ signature and the lowest innate anti-PD-1 resistance signature among subtypes. TP53-wt revealed a converse enrichment trend for these immunologic signatures. Survival analyses using the Cancer Genome Atlas Uterine Corpus Endometrial Carcinoma (TCGA-UCEC) dataset identified associations of CCR5 (hazard ratio (HR) = 0.71, p = 0.035), TNFRSF14 (HR = 0.58, p = 0.028), and IL-10 (HR = 2.5, p = 0.012) with overall survival within MSI-H, TP53-mutant, and TP53-wt subtype, respectively. A TIME comparison between the sub-stratified subgroups of our cohort revealed upregulated tumor infiltration of immune cells in the low-risk subgroups. Our study demonstrates that targeted DNA sequencing is an effective one-stop strategy to classify EC. Immunomodulatory genes may serve as prognostic markers within subtypes.


Introduction
Uterine corpus cancer is the second most common gynecological cancer worldwide [1] and its incidence is rapidly rising in China [2]. Endometrial cancer (EC) comprises most uterine cancer cases. As a highly heterogeneous entity, EC can be divided into various histologic subtypes: mainly endometrioid endometrial cancer (EEC), serous endometrial cancer (SEC), and clear cell endometrial cancer (CCEC). Currently, pathologic evaluation is routinely used in clinical practice to predict prognosis and guide the management of EC patients.
With the expanded understanding of the EC genome, extensive investigations have demonstrated the potential clinical utility of molecular subtyping EC for risk stratification and selecting optimal treatment. Using array and sequencing-based technologies, the Cancer Genome Atlas (TCGA) research network first divided EECs and SECs into four molecular subtypes with prognostic significance: the POLE/ultramutated subtype exhibited the most favorable prognosis; the microsatellite instability (MSI)/hypermutated and copy-number low subtypes showed intermediate prognosis; the copy-number high subtype, also characterized by a low mutation burden and a high TP53 mutation frequency, had poor prognosis [3]. Afterward, to broaden the utility of TCGA classification, other research groups used immunohistochemical (IHC) assessment of mismatch repair (MMR) proteins and p53 as surrogates for testing MSI and somatic copy-number alterations (SCNA). The Translational Research in Post-Operative Radiation Therapy in EC (TransPORTEC) system defines four subtypes: p53 abnormal, MSI, POLE-mutant, or no specific molecular profile (NSMP) [4,5]. The ProMisE (Proactive Molecular Risk Classifier for Endometrial Cancer) system classifies patients in the light of testing for aberrations in the sequence of MMR deficiency (dMMR), POLE mutation, and p53 abnormality [6]. Despite showing promising values in risk stratification, these molecular classifiers integrate results from multiple testing platforms that are labor-intensive. A highly professional and multi-disciplinary team is also required to interpret the results. These drawbacks will largely confine the clinical application of these classification schemes. Developing an easy-handling and clinically feasible approach for molecular subtyping of EC remains an unmet need.
Although incorporating molecular subtypes into risk stratification can limit the current under-and overtreatment of EC patients, some molecular subtypes still lack a univocal prognostic significance, especially the NSMP group, the prognosis of which is largely affected by clinicopathological features of tumors [4,5]. The intra-subtype heterogeneity suggests that further stratifying subtypes based on prognostic factors may improve the current molecular classification [7].
Mounting evidence reveals the role of the immune system in the progression and prognosis of various tumors, including EC. The favorable prognosis of POLE-mutant and dMMR subtypes has been suggested to correlate with increased density of CD8+ tumorinfiltrating lymphocytes (TILs) and higher PD-1expression [8,9]. van Gool et al. showed that the upregulated cytotoxic T cells in POLE-mutant tumors were accompanied by enriched tumor-infiltrating T-cell gene expression signature and increased expression of T-cell cytotoxic differentiation and effector markers [10]. Talhouk et al. characterized the immune microenvironment across the four molecular subtypes of EC and found profound variation in immune response across and within subtypes [11]. Collectively, these results suggest the potential for using immunologic signatures to sub-stratify patients within subtypes.
Herein, we classified women with EC using a one-stop strategy based on targeted DNA sequencing and interrogated immunologic signatures across subtypes, including transcriptome of immunomodulatory genes and tumor immune microenvironment (TIME) landscape. We also investigated the prognostic value of immunomodulatory genes within individual subtypes using the TCGA cohort.

Characteristics of Patients
A total of 50 women with EC were included in the study, with a median age of 57 years. Most patients presented with stage I or II (62%), endometrioid (86%), and low-grade (G1 or G2, 70%) tumors. Three and four patients had serous and clear cell histology, respectively ( Table 1). The cohort had a postoperative follow-up of 10 to 14 months. Up to the last follow-up date, the progression-free survival (PFS) rate was 94%. Three (6%) out of fifty patients developed recurrent or progressive disease.
The TP53-mutant subtype showed lower estrogen receptor (ER) and progesterone receptor (RP) positive rates than other subtypes (p = 0.006, p = 0.002, Table 1). A higher frequency (100%) of lymphovascular space invasion was observed in the POLE-mutant subtype (p = 0.048). TP53-wt tumors displayed lower Ki67 expression than other subtypes (p = 0.001). We also performed IHC for MMR and p53 proteins and classified patients according to the ProMisE system. Figure 1B illustrates the comparison of classification between the two approaches. Of the 50 patients, only 5 had discordant subtypes, resulting in a concordance of 90%. Specifically, patients p06 and p07 were identified with an MSI-H status by NGS but with a proficient MMR (pMMR) and p53-wt status by IHC, therefore were classified into the subtype of MSI-H by NGS and p53-wt by ProMisE. Patients p21 and p24 were determined as microsatellite stable (MSS) and TP53-mutant by NGS but deficient MMR (MMR-D) by IHC, resulting in the discordant subtyping between NGS (TP53-mutant) and ProMisE (MMR-D). Patient 41 was classified as TP53-wt by NGS and p53-abn by ProMisE.

Differential Immunologic Signatures across Four Subtypes
Next, we interrogated the transcriptome profile of 83 immunologic genes among four EC molecular subtypes ( Figure 3A). Differential expression levels were observed across subtypes for numerous genes related to T-effector and interferon-gamma (IFN-γ) gene signature and innate anti-PD-1 resistance and CD8+ T cell exhaustion (IPRES_eCD8T) signature (Table S1). When comparing signature enrichment based on single-sample Gene Set Variation Analysis (ssGSVA) scores, we found POLE-mutant tumors showed the highest enrichment for the signature of T-effector and IFN-γ, and the MSI-H subtype ranked second, followed by TP53-wt. TP53-mutant tumors did not show significantly differential enrichment for this signature vs. other subtypes ( Figure 3B). Conversely, the enrichment score for the IPRES_eCD8T signature was lower in the POLE-mutant vs. TP53-wt tumors (p < 0.01), and the latter showed a higher score than the TP53-mutant subtype (p = 0.025) ( Figure 3C). We also observed higher enrichment for the epithelial-mesenchymal transition (EMT) signature in TP53-wt tumors than in TP53-mutant (p < 0.01) and MSI-H subtypes (p = 0.014) ( Figure 3D). The transforming growth factor-Beta (TGF-Beta) signature was not enriched differentially across four subtypes ( Figure 3E). We also investigated the TIME landscape among different subtypes by assessing the tumor-infiltrating immune cells. A differential rate of CD8+ TILs was observed across subtypes in the tumor epithelial compartment (p = 0.029, Table S2). Specifically, POLEmutant and TP53-mutant tumors showed higher rates of epithelial CD8+TILs than TP53-wt tumors ( Figure S2A). In the stromal compartment, POLE-mutant (p < 0.01, p < 0.01) and MSI-H (p < 0.05, p < 0.01) subtypes expressed both higher densities and rates of PD-L1+ marker than TP53-wt tumors, and a higher rate of CD8+TILs was observed in POLE-mutant (p < 0.05) and TP53-mutant (p < 0.01) tumors than in TP53-wt subtype ( Figure S2B). Other immune markers were generally expressed comparably across different molecular subtypes (Table S2).

Prognostic Value of Immunomodulatory Genes
Despite differential expression profiles of immunomodulatory genes across subtypes, we also observed heterogeneous expression patterns among tumors within the same subtype ( Figure 3A), suggesting the potential of further stratifying patients using these genes. Therefore, we retrieved mRNA expression data of the 83 genes covered in the present study and clinical information from 232 EC patients, whose molecular subtypes had been released [3], from the TCGA Uterine Corpus Endometrial Carcinoma-(TCGA-UCEC) database. First, we analyzed the association of individual gene expression with prognosis and identified 24 genes significantly associated with overall survival (OS) ( Table S3, p < 0.05). Subsequently, we assessed the prognostic value of these candidate genes within each molecular subtype, except for POLE-mutant tumors, which invariably exhibit the most favorable prognosis regardless of clinicopathological features [7]. Using the median expression value of CCR5, we stratified patients with MSI-H tumors into two subgroups with differential OS (hazard ratio (HR) = 0.71, p = 0.035) ( Figure 4A). Conversely, a higher expression of IL10 yielded a shorter median OS within the TP53-wt subtype (HR = 2.5, p = 0.012, Figure 4B). Within TP53-mutant subtype, we found higher expression levels of CCR5 (p = 0.003), CD2 (p = 0.037), CXCL11 (p = 0.033), LAG3 (p = 0.045), TBX21 (p = 0.014), TNFRSF14 (p = 0.004), and TIGIT (p = 0.006) were associated with more favorable survivals ( Figure 4C). Multivariate analysis confirmed TNFRSF14 as an independent prognostic factor in patients with TP53-mutant EC (HR = 0.58 (0.36−0.94), p = 0.028).
Next, we sub-stratified the patients of our cohort using the same cut-offs of the abovementioned prognostic genes and compared the TIME landscape between subgroups within each molecular subtype (Table S4). The CCR5-high subset of MSI-H tumors revealed higher positive rates of intraepithelial CD3+ and CD8+ TILs compared with the CCR5-low subset ( Figure 5A). The TNFRSF14-high subset of TP53-mutant tumors expressed higher positive rates of intraepithelial PD-L1+ markers and CD8+ TILs than the TNFRSF14-high subset ( Figure 5E). TIME signatures were generally comparable in the stromal compartment between subsets of MSI-H and TP53-mutant tumors ( Figure 5B,F). Moreover, we did not observe significantly differential TIME signatures in the epithelial compartment between the two subsets of TP53-wt tumors sub-stratified by IL10 expression (Figure 5C). While the positive rate of CD163−CD68+ cells (putative M1 macrophages) in the stromal area tended to be lower in the IL10-high subset (p = 0.078, Table S4, Figure 5D)).  Notably, one patient with stage IIIC1, G1, and MSI-H EC developed disease progression during the postoperative treatment with radiotherapy, chemotherapy, and PD-1 inhibitor. The patient was sub-stratified as CCR5-low, suggestive of an inferior prognosis. In the TP53-wt group, one patient with stage IIIC1 G1 EC refused adjuvant therapy, and systemic recurrence was noted six months after surgery. The patient was also sub-stratified into a high-risk subset featured by IL10-high. In the TP53-mutant group, one patient with high-grade serous carcinoma staged as IVB experienced intraperitoneal recurrence two months after finishing the standard adjuvant therapy. The patient was characterized as TNFRF14-low, indicating a poor prognosis.

Discussion
In this study, we performed molecular subtyping of 50 EC patients using targeted DNA sequencing and discovered differential immunologic gene signatures and TIME across varied subtypes. In addition, we also identified CCR5, TNFRSF14, and IL-10 that are associated with OS within MSI-H, TP53-mutant, and TP53-wt subtypes, respectively.
Our study classified 43 EECs, 3 SECs, and 4 CCECs into four molecular subtypes through a one-stop NGS-based strategy. Compared with TransPORTEC and ProMisE systems, our approach, which can be conducted largely in an automatic manner, is more labor-saving. It also relies less on experienced pathologists, which may introduce interobserver variability. In our study, four patients show discordant results between IHC and NGS results for MSI. Two MSI-H/pMMR cases can be caused by inactive mutant proteins that remain detectable by IHC. Due to functional redundancy, cases with isolated loss of MMR protein may still have MMR function, which may explain the MSS/dMMR results in two patients [12,13]. In these two scenarios, NGS outperforms IHC since it reflects the actual status of MMR function. Moreover, because it was found in 92% of copy number-high tumors, the TP53 mutation is considered a surrogate for copy number-high subtyping [3]. Given the discordance of~8% between mutation and protein status of TP53, the routinely used IHC would misclassify~15% copy number-high tumors into p53-wt/NSMP subtype [14,15]. In this regard, our NGS strategy would reduce such misclassifications. However, it should be noted that the current cost of NGS may prevent the implementation of our strategy. However, with the rapid development of high throughput sequencing, the cost per sample will ultimately reduce to an affordable range. In addition to molecular subtyping, NGS can provide information that may guide the therapeutic decision. This advantage will promote the utility of the NGS strategy in the clinic.
We present the first study comprehensively investigating the immunomodulatory transcriptome signatures across four molecular subtypes of EC. Our study showed an upregulation of T-effector and IFN-γ signature in POLE-mutant and MSI-H tumors and downregulation of IPRES signature in POLE-mutant tumors. The gene signature of Teffector and IFN-γ, indicative of a T-cell-activated tumor environment, has been associated with response to PD-1 blockade [16]. Contrastingly, the IPRES signature defines a set of genes related to immune-suppressive cytokines, EMT transcription factors, and proangiogenic factors and has been associated with innate resistance to PD-1 blockade across multiple cancers [17]. Moreover, the TP53-wt subtype expressed a higher level EMT signature. Despite being a major inductor of EMT, TGF-β related gene expression was comparable among different subtypes. This discordance may be due to the relatively limited sample size of our cohort, given that the TP53-wt tumors showed a numerically higher median level of TGF-β signature than TP53-mutant and MSI-H tumors ( Figure 3E). Another possible explanation is that the EMT process in the TP53-wt tumors can be induced by activation of other signaling pathways such as the Wnt/β-catenin pathway [18].
Our study also depicts differential TIME landscapes across four subtypes of EC. Higher rates and densities of PD-L1+ cells were found in stromal areas of POLE-mutant and MSI-H tumors, and higher rates of CD8+TILs were seen in both tumoral and stromal areas of POLEmutant tumors. We also found lower rates of both tumoral and stromal CD8+TILs present in TP53-wt than in TP53-mutant tumors. Consistently, studies have described higher numbers of CD8+ TILs in POLE-mutant and MSI-H ECs than in MSS tumors [8,9]. Analysis of the TCGA EC data set revealed that POLE-mutant tumors express higher transcripts of PD-L1, PD-L2, CD8A, and PD-1 [19]. Talhouk

et al. found elevated infiltration by both T cells (CD3+CD8+ TILs and CD3−CD8+ TILs) and B-lineage cells (CD138−CD79a+ and
CD138+CD79a+) in POLE-mutant and MMRd tumors than in p53-abn and -wt subtypes. The authors also demonstrated a lower level of CD8+ TIL activation and/or exhaustion in p53wt tumors than the other three subtypes [11].
MSI-H and TP53-wt subtypes, accounting for approximately 70% of ECs, display intermediate and indistinguishable prognoses [4,5], highlighting the need for the sub-risk stratification of patients within these two subtypes. Intriguingly, we identified several immunological genes with potential prognostic value within molecular subtypes of EC. High CCR5 expression was a predictor for longer OS within MSI-H tumors. CCR5 gene encodes C-C chemokine receptor type 5, a seven-transmembrane G protein-coupled receptor that binds to various ligands. The CCL5/CCR5 axis promotes tumor progression through a spectrum of mechanisms, such as enhancing invasion and metastasis of tumors, reducing tumor cell resistance to drugs, and recruiting immune and stromal cells [20]. CCR5 is selectively overexpressed in breast cancer cells and promotes tumor metastases, which is associated with poor prognosis [21]. In colorectal cancer cells, the blockage of CCR5 suppresses cell proliferation, migration, and clonogenic ability [22]. On the other hand, as a double-edged sword, the CCL5/CCR5 axis also promotes antitumor immunity by recruiting tumor-infiltrated-T cells and dendritic cells [23,24]. Our study demonstrated CCR5 as a predictor of favorable prognosis in MSI-H EC for the first time. We also found that the upregulation of CCR5 was accompanied by elevated intraepithelial CD3+ and CD8+ TILs in MSI-H EC. Our observation is consistent with the previous finding that CCR5 expression is positively associated with tumor immune cell infiltration [25]. Moreover, Talhouk et al. described a positive association of CD8+CD3+and CD3+CD8− TILs with disease-specific survival(DSS) within MSI-H EC [11]. Collectively, all these observations support the prognostic significance of CCR5 expression within MSI-H EC.
Our study identified IL-10 as a potential adverse prognostic factor within TP53-wt tumors. IL-10 gene encodes for interleukin 10, a cytokine that serves as a pleiotropic immunomodulatory factor, predominantly with immunosuppressive effects [26,27]. The increased serum level of interleukin 10 has been described as a predictor of unfavorable outcomes in patients with several cancers (colorectal cancer, renal cell carcinoma, etc.) [28][29][30]. Our study is the first one demonstrating the adverse prognostic effect of IL-10 gene expression in TP53-wt EC. Our results also revealed that the low-IL10 subgroup (predicted to have a better prognosis) tended towards a TIME landscape of upregulated stromal M1 macrophages compared with the high-IL10 subgroup. Interestingly, Talhouk et al. showed an association of stromal B cells with increased disease-specific survival within the TP53-wt EC [11]. Activated M1 macrophages release pro-inflammatory factors and suppress tumor growth. High infiltration of M1 macrophages has been described in EC patients with a better prognosis [31].
Our study also identified TNFRSF14 as an independent prognostic factor within the TP53-mutant subtype, with higher expression associated with longer OS. TNFRSF14 encodes for a receptor (HVEM) expressed on T cells, B cells, monocytes, and immature dendritic cells. It functions in lymphocyte activation and regulates the antitumor immune response [32]. TNFRSF14 induces apoptosis and suppresses proliferation, thus playing a tumor-suppressive role in bladder cancer, which may confer a better prognosis for patients with upregulated TNFRSF14 expression [33]. In vitro evidence revealed that the inactivation of TNFRSF14 promotes the migration and invasiveness of EC cells by inducing the EMT [34]. Our results, for the first time, demonstrate the prognostic significance of TNFRSF14 within the TP53-mutant subtype. The prolonged survival of the subgroup expressing higher TNFRSF14 is likely attributed to the increased intraepithelial CD8+ TILs ( Figure 5E), which has been previously associated with a better prognosis in EC patients [35]. However, we also observed more PD-L1+ cells in the high TNFRSF14 subgroup. The association of PD-L1 expression with prognosis remains elusive in EC patients since relevant studies showed controversial results [36,37]. Our study suggests overexpression of PD-L1 may impact patient outcomes in a heterogeneous manner across varied subtypes of EC, which merits further exploration.
The major limitation of our study is that the follow-up time is not long enough to obtain mature survival data; thus, the prognostic significance of genes identified from the TCGA dataset cannot be validated in our cohort. The differential TIME landscapes between subgroups of our cohort stratified by transcriptome signatures might suggest their prognostic impact. Furthermore, the study enrolled a relatively small number of patients. After molecular classification, the number of patients in each subtype was limited. The finding of our study definitely needs further validation in a larger external cohort with a longer flow-up time.
In conclusion, our study demonstrates that targeted DNA sequencing is an effective one-stop strategy to classify EC patients. Moreover, we comprehensively interrogated the immunomodulatory transcriptome signatures and TIME landscape across four subtypes. Our study identified several immunomodulatory genes that may serve as prognostic predictors to sub-stratify patients within subtypes. Our results reveal the potential of refining the risk stratification in EC using the immunomodulatory signature, which may better guide the clinical management and decision-making for patients with EC. As illustrated in Figure 6, DNA sequencing was performed with 50 formalin-fixed, paraffin-embedded (FFPE) tumor samples for molecular subtyping of patients according to POLE, TP53, and MSI status. MLH1, PMS2, MSH2, MSH6, and p53 were evaluated using IHC staining. TP53 and MSI statuses determined from DNA sequencing were compared with corresponding IHC results. RNA sequencing and multiplex immunofluorescence assay were conducted in parallel with the same 50 samples to compare the expression profile of 83 immune-related genes and TIME landscape among different molecular subtypes. Gene expression (https://portal.gdc.cancer.gov/projects/TCGA-UCEC, accessed on 5 June 2021) and survival data (https://xenabrowser.net/datapages/, accessed on 5 June 2021) of TCGA-UCEC were retrieved to screen for prognostic genes. The prognostic value of selected candidates was further investigated within each subtype. Subsequently, patients of each subtype from our cohort were sub-stratified according to the expression of the identified prognostic gene, and the TIME landscape was compared between subgroups. All procedures performed in studies involving human participants were in accordance with the 1964 Helsinki declaration and its later amendments (https://www.wma.net/policies-post/wma-declaration-of-helsinki-ethicalprinciples-for-medical-research-involving-human-subjects/, accessed on 5 June 2021). The study was approved by the Institutional Ethics Committee at Peking Union Medical College Hospital (No. JS-2884). All patients provided written informed consent.

DNA Sequencing and Molecular Classification of Patients
DNA was extracted from FFPE tumor tissues using the QIAamp DNA FFPE tissue kit (Qiagen, Hilden, Germany), following the manufacturer's instructions. Capture-based targeted sequencing was performed using a 520-gene panel (OncoScreen ® Plus, Burning Rock Biotech, Guangzhou, China) as previously described [38][39][40]. Data analyses, including variants calling and interpretation, copy number variation, TMB estimation, and MSI status assessment, were carried out using standardized pipelines based on the methods described previously [40].
Patients were subtyped based on the genomic signatures identified by DNA-sequencing in the order of POLE hotspot mutation (P286R, V411L, S297F, A456P, or S459F), MSI-H, and TP53 mutation. Patients without any of the three signatures were determined as TP53-wt.

RNA Sequencing and Gene Expression Profiling
RNA was extracted from FFPE samples using an All prep DNA/RNA FFPE kit (Qiagen, Hilden, Germany). The quantity and quality of extracted RNA were measured using Qubit RNA HS assay (Thermo Fisher Scientific, Waltham, MA, USA) and LabChip GXII touch 24 (Perkin Elmer, Waltham, MA, USA), respectively. RNA sample with a concentration ≥ 4ng/uL and DV200 ≥ 40% was considered eligible for subsequent library construction. After strand-specific cDNA synthesis, dA-tailing, unique molecular identifier adaptor ligation, and PCR amplification, the products were hybridized with the capture probe of the OncoRNA panel (Burning Rock Biotech, Guangzhou, China). The panel consists of 83 immunomodulatory genes (Table S1). The prepared libraries were sequenced, and the sequencing data were processed as previously described [41]. Gene-level expression values were computed as transcripts per million (TPM) and normalized to Z-scores before clustering. A heatmap of gene expression was drawn using Pheatmap (1.0.12). Genes were categorized into different gene sets designated as signatures of T-effector and IFN-γ, IPRES_eCD8T, EMT, and TGF-Beta (Table S1), according to previous publications [17,42,43].
The ssGSVA program [44] was used to calculate a single-sample gene set enrichment score for each signature. A total of 50 samples had been finally analyzed.

Multiplex Immunofluorescence Assay
FFPE blocks were serially sliced into sections of 5 µm and subjected to multiplex immunofluorescence staining using the PANO 7-plex IHC kit (Panovue, Beijing, China) according to the manufacturer's instruction. Two panels were analyzed in this study and stained with different combinations of primary antibodies: one was stained with antibodies for PD-L1, PD-1, CD3, CD8, and Pan-CK and the other was stained with antibodies for CD68, CD163, CD56, and Pan-CK. Subsequently, slides were incubated with horseradish peroxidase (HRP)-conjugated secondary antibody followed by tyramide signal amplification using TSA Fluorescence Kits (Panovue, Beijing, China).
The Mantra System (PerkinElmer, Waltham, MA, USA) was used to scan the slides. The fluorescent spectra were set at 20-nm wavelength intervals from 420 to 720 nm. A spectral library was established using the extracted images and used for multispectral unmixing by inForm image analysis software (PerkinElmer, Waltham, MA, USA). The present study reported the densities and positive rates of PD1+, PD-L1+, CD3+ (designating T cells), CD8+ (designating cytotoxic T cells), CD56+ (designating for NK cells), CD163−CD68+ (designating M1 macrophages), CD163+CD68+ (designating M2 macrophages) markers on cells in both the epithelial and stromal compartments of tumors.

Statistical Analyses
All statistical analyses were performed using R software (version 3.5.3). Differences among groups were compared using the Chi-square test or Fisher's exact test for categorical variables and the Kruskal-Wallis test for continuous variables. The Wilcoxon test was used to compare to differences between groups. Multivariate Cox regression analysis was performed to adjust for confounders. Kaplan -Meier analysis was used to estimate survival, and a log-rank test was used to determine the differences in the multiple survival metrics between groups. Statistical significance was defined as p < 0.05. Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The datasets generated and/or analysed during the current study are available in the NGDC repository [https://ngdc.cncb.ac.cn/, OMIX002704, OMIX002705].

Acknowledgments:
The authors would like to thank Wang Wenze and Cuijie from Department of Pathology, Peking Union Medical College Hospital for conducting part of the histopathological experiments and confirming the results. The authors thank Lin Shao, Danhua Wang, Lan Su, and Qiang Lu from Burning Rock Biotech for their help in data analyses, language editing and manuscript formatting.

Conflicts of Interest:
The authors declare no conflict of interests.