Dissecting the tumour immune microenvironment in merkel cell carcinoma based on a machine learning framework

Abstract

Merkel cell carcinoma (MCC) is a primary cutaneous neoplasm of neuroendocrine carcinoma of the skin, which is characterized by molecular heterogeneity with diverse tumour microenvironment (TME). However, we are still lack knowledge of the cellular states and ecosystems in MCC. Here, we systematically identified and characterized the landscape of cellular states and ecotypes in MCC based on a machine learning framework. We obtained 30 distinct cellular states from 9 immune cell types and investigated the B cell, CD8 T cell, fibroblast, and monocytes/macrophage cellular states in detail. The functional profiling of cellular states were investigated and found the genes highly expressed in cellular states were significantly enriched in immune- and cancer hallmark-related pathways. In addition, four ecotypes were further identified which were with different patient compositions. Transcriptional regulation analysis revealed the critical transcription factors (i.e. E2F1, E2F3 and E2F7), which play important roles in regulating the TME of MCC. In summary, the findings of this study may provide rich knowledge to understand the intrinsic subtypes of MCCs and the pathways involved in distinct subtype oncogenesis, and will further advance the knowledge in developing a specific therapeutic strategy for these MCC subtypes.

Keywords:

• Tumor microenvironment
• machine learning
• cellular states
• merkel cell carcinoma
• transcriptional regulation

Introduction

Merkel cell carcinoma (MCC) is a rare but highly aggressive skin cancer, which has neuroendocrine features and a high local recurrence rate and the potential for distant metastasis [1]. Indeed, MCC is one of the high mortality rate of cancers, with more than one-third of MCC patients dying and is twice as deadly as malignant melanoma [2]. It was discovered that a double-stranded DNA polyomavirus was integrated into the genomes of eight tumour cells in MCCs, suggesting that polyomaviruses are oncogenic [3]. Moreover, it has been shown that the vast majority of MCCs are related to the clonal integration of the Merkel cell polyomavirus (MCPyV), and only a small number of MCCs are due to a large number of mutations caused by prolonged sun exposure, which are known as MCPyV-MCC [3], suggesting that Merkel cell polyomavirus infection was the main risky factor for MCC. Additionally, MCPyV + MCC patients had a three-fold higher survival rate than those with MCPyV- MCC, MCPyV can be used as prognostic factor [4]. Moreover, the tumour immune microenvironment (TME) plays an important role in various types of cancer [5,6], but the in-depth study of TME in MCC is still seriously lacking.

The TME is composed of many different cell types, including B cells, fibroblasts, T cells, and dendritic cells [7]. Thus, it is essential to analyze the different cell types in TME of MCC. Recent studies have demonstrated that MCPyV + MCC has higher immune infiltration [8,9]. Considering the TME of MCC, previous studies mainly focused on the expression of PD-L1 in tumour-infiltrating lymphocytes [8,10], and have shown that significantly different in PD-L1 expression in tumour cells due to the presence of MCPyV [8]. In addition, characterizing the TME of MCC was performed by comparing the expression of important genes. For example, CD200 can promote the development and progression of a variety of cancers [11–13]. In addition, CD200 was frequently and highly expressed in MCC, suggesting that CD200 may be a novel immune target independent of T-cell immune checkpoints [14]. Several lines of evidences suggested several other genes were significantly differently expressed between MCPyV + MCC and MCPyV- MCC [15], including TP53 and RB1. The expressions of these genes were associated with clinical survival of cancer patients [16–18].

Moreover, recent cancer transcriptome analysis has revealed different immune subtypes of cancer [19,20], which extends the conventional cancer classification based on anatomical site of original [21]. Transcriptome-based classifications of tumour subtypes have informed a number of characteristics of clinical patients, including the response to immunotherapy. But it usually oversimplifies the cell types in TME and not considered the cellular states and ecosystems of tumours [22]. Recent developments of single-cell RNA sequencing (scRNA-seq) have enabled comprehensively surveys of TMEs of various tumour types [23]. However, there is still lack of scRNA-seq data for MCC. Recently, Luca, Steen et al. develop the EcoTyper software to deconvolve tumour-microenvironment interactions from high volume bulk transcriptomics data [24]. This tool has been successfully applied in various types of cancer, such as diffuse large B cell lymphoma (DLBCL) [25]. However, the landscape of cellular states and ecosystems in MCC is still unknown.

Here, we comprehensively characterized the cellular states and ecosystems in MCC by EcoTyper, which is a machine learning framework. We resolved the immunophenotype of MCC based on 12 cell types. We obtained the landscape of cellular states and ecotypes, which were enriched the existing immunotypes. In addition, marker genes of cellular states were analyzed to understand the biological pathways affected by each cell state and we revealed the biomarkers of MCC. The potential transcription factors (TFs) that regulating the cellular states were also identified. Together, our comprehensive transcriptome integrative analyses provide a comprehensive atlas of TMEs in MCC and functional characterization of TMEs greatly helps identify biomarkers of MCC.

Materials and methods

Transcriptomes of MCC

We performed a systematic analysis of public gene expression data of Merkel cell carcinoma (MCC). The transcriptomes of MCC patients were downloaded from gene expression omnibus (GEO) under the accession number GSE39612 [16]. The datasets were obtained from the Affymetrix® GPL570 platform (Human Genome U133 Plus 2.0 Array). We downloaded the raw CEL file from GEO, and processed via the robust multiarray averaging (RMA) algorithm implemented in the ‘affy’ package in R program. All probes were mapped to gene symbols by the annotation in ‘hgu133plus2.db’, and the probes that mapped to multiple gene symbols were removed. Finally, there were 20,862 genes remained for further analysis and 74 tumour patients of MCC were considered in this study.

Cell-type specific gene signatures

To obtain the cell-type specific gene expression profiles of MCC patients, we used CIBERSORTx [26] to estimate gene expression profiles and assess the proportion of major cell populations. To conduct a comprehensive analysis of MCC, we applied the signature matrix published in the EcoTyper [22], which contains 12 major cell types. This signature matrix is a combination of two previously published signature matrices: LM22, which includes 22 human immune cell types, and TR4, which includes epithelial, endothelial, immune and fibroblast populations.

To obtain the signature matrix suitable for cell classification, the cell types in the LM22 were mapped into 9 major lineages: B cells, plasma cells (PC), CD8 T cells, CD4 T cells, natural killer (NK) cells, monocytes/macrophages, dendritic cells, mast cells, polymorphonuclear neutrophils (PMNs). Based on the transcriptomes of MCC, we obtained fractional abundance estimates for 12 major cell types. In each sample, the summarized abundances 12 cell types were normalized to 1. Finally, CIBERSORTx was employed to determine gene expression profiles by default parameters.

Discovery of cellular states in MCC

After obtained the cell-type specific gene expression profiles, we employed the machine learning method, EcoTyper, to classify clusters for each cell type. Based on non-negative matrix factorization (NMF) combined with specific heuristics, EcoTyper calculated the generated cell-type specific gene expression profiles to identify and quantify cellular states for each cell type. In the process of calculation, EcoTyper determined the number of cell states for each cell type based on cophenetic coefficient and the cluster number was default ranged from 2 to 20. We selected the cluster number closest to a cophenetic coefficient of 0.99 as the most suitable cluster number in each cell types. Cellular states whose marker genes less than 10 and that likely to be false positives according to false-positive index (AFI) were removed in our analysis.

Profiling of ecosystems in MCC

Equally important, we employed EcoTyper [22] to identify the robust cellular state co-occurrence network terms ecotypes. EcoTyper uses a community identification approach that the pairing co-associations between cell states were maximized while mutual repulsive forces were minimized. In brief, a denote binary matrix with cellular states as rows and samples as columns, if cell state i is assigned to sample j, A^i,j is set 1, otherwise is 0. Next, the Jaccard matrix of paired combinations of all cell states in A was calculated. According to the number of MCC cellular states, a 30*30 Jaccard matrix was generated. After the Jaccard matrix was generated, the hypergeometric test was carried out on pairing cell states (such as s and p), entries in J^s,p are set 1 if p is greater than 0.01 and 0 otherwise.

Next, in order to identify ecotypes, ‘hclust’ function in R package stats was used to conduct unsupervised clustering of Jaccard matrix. The optimal number of clustering was determined according to silhouette width maximization. Finally, two ecotypes were obtained in MCC. Cytoscape [27] was then used to generate network diagrams of MCC ecotypes.

Functional analysis of cellular states

To examine functional landscape of cellular states and ecotypes, we performed functional enrichment analysis for marker genes of cellular states and ecotypes obtained from EcoTyper. First, we used ‘clusterProfiler’ [28] package to analyze biological processes in Gene ontology (GO). Only the GO terms with genes ranges from 15 to 500 were considered, the significance thresholds of biological processes were p < 0.01 and adjust p < 0.05. And for each cell state, we show the visualization of top ten biological functions the highest gene counts. The GO terms of ecotypes were clustered based on simplify Enrichment R package [29]. The similarities among GO terms were calculated by the ‘GO_similarity’ function and the cluster results were visualized by ‘simplifyGO’ function.

To analyze the effects of each cellular state on the cancer hallmark pathways and the immune-related pathways, we downloaded the cancer hallmark pathways from MsigDB [30] and the immune-related pathways were downloaded from ImmPort project [31]. We used hypergeometric test to analyze these pathways and we only considered pathways with adjust p < 0.05 as significant. Moreover, we applied the Gene Set Enrichment Analysis (GSEA) to further analyze genes highly expressed in each cellular state [32].

Transcriptional regulation analysis

To determine the transcription factors (TFs) that potentially regulated the highly expressed genes in ecotypes, we performed the transcriptional regulation analysis by the RcisTarget. First, we used Wilcoxon’s rank sum test to obtain the highly expressed genes in each ecotype. The genes with adjusted p-values less than 0.05 were selected as the highly expressed genes in each ecotype, and applied these genes to RcisTarget. We applied the TF binding sites upstream 500 bp of transcription start sites and ‘cistarget’ function was used to analyze with the parameters (nesThreshol = 3, AUC = 3*sd(auc) + mean(auc)). The enriched motifs were ranked based on the normalized enrichment score (NES) and top five were visualized for each ecotype in MCC.

Results

High-throughput cellular states discovery in MCC

We applied EcoTyper to identify cellular states from bulk transcriptome profiles of MCC patients. The default first step of EcoTyper is to estimate gene expression profiles and assess the proportions of major cell populations with CIBERTSORTx. To compare the cellular heterogeneity from 74 MCC patients, the abundance of cell type was used to estimate cell-type specific gene expression profiles of 12 cell types. Next it identified cellular states transcriptionally defined and confirmed robust cellular state co-occurrence that defined the cellular communities (Figure 1A).

Figure 1. The machine learning framework for high-throughput characterization of cellular states and ecosystems in MCC. (A) Overview of the application of EcoTyper to cellular states and ecosystems profiling in MCC. (B) UMAP projection of cellular states across MCC tumour patients. (C) The proportions of patients within different cellular states in each cell type. (D) The sankey diagram showing the relationship of epithelial cellular states, tissue, and merkel cell polyomavirus status.

Based on the computational analysis, we found that 74 MCC patients were divided into 9 cellular types, including B cells, CD8 T cells, Dendritic cells, Endothelial cells, Epithelial cells, Fibroblasts, Mon/Mac, PCs, and PMNs. Next, we obtained 30 distinct cellular states from these 9 cell types (Figure 1B). The numbers of cellular states were clearly different, ranging from 2 to 10 in each cell type. Moreover, two or three cellular states were observed in the majority of cell types, but there were ten clusters in Endothelial cells (Figure 1B).

To clearly show the distribution of samples for cellular states in each cell type, we calculated the proportions of patients. We found that the proportions of patients were significantly different in each cellular state (Figure 1C). For example, patients were mainly divided into S01 in the Fibroblasts and PCs, whereas S01 had fewer patients than the other states in Dendritic cells. Recent studies have shown that MCC presents a trilinear different characteristics, concurrenting the expression of B-lymphoid lineage, epithelial and neuroendocrine markers [33]. In most studies it has been pointed out that epithelial stem cells are the most likely cellular origin of MCC [34–36]. For epithelial cells, we obtained two celluar states (Figure 1B). By analysis of the tissue and virus infection status of MCC patients in each cellular state, we found that there were more primary MCC patients in S01, and the majority of patients in S01 were not infected with MCPyV (Figure 1D). Together, these results provided a landscape of cellular states in MCC, which will be a valuable resource for understanding the TME of MCC patients.

B cell and CD8 T cell states in MCC

To further analyze the cellular state profile of MCC, we started with specific cell types. Recent finding indicated that cellular immunity plays a key role in MCC [37]. Meanwhile, increasing evidences have shown that MCCs express the marker genes of early B cell [18,38]. B cells were hypothesized to be the origin cells of MCC by the consistency of early B cell antigen expression in MCC [39]. Therefore, it is vital to delineate the subtypes defined by transcription of these immune cells. For B cells in MCC, we identified two cellular states by EcoTyper (Figure 2A), whereas we obtained three cellular states in CD8 T cells (Figure 2B). We presented the samples information of CD8 T cells and B cells in detail, and found that most cases in S01 of both cell types are primary MCC, and most of them are MCPyV + MCC, this was consistent with the conclusion that MCPyV + MCC patients have better survival (Figure 2C–D).

Figure 2. The cellular states landscape of B cell and CD8 T cell in MCC. (A-B) Heatmap depicting cellular states of B cells and CD8 T cells identified from MCC bulk transcriptomes. Patient samples are organized by the most prevalent cell states and genes used for identification of the cell states are shown. (A) for two B cell states. (B) for three CD8 T cell states. (C-D) a sankey diagram showed the relationship of cellular states, tissue, and merkel cell polyomavirus status. (C) for B cell states. (D) for CD8 T cell states. (E) Boxplot showing the expressions of WDFY4, CD19, CD38 and MS4A1 in patients within B cell states. (F) Boxplot showing the expressions of IL16, IL2RB, IL18BP and MPEG1 in patients within CD8 T cell states. (G) Barplots showing the enriched functional pathways for genes highly expressed in B cell S01 and CD8 T cell S01.

In addition, we found that genes involved in immune response were specifically highly expressed in the corresponding cellular states (Figure 2A). For example, WDFY4, CD38, CD19, MS4A1/CD20 were highly expressed in S01 state of B cells (Figure 2E). WDFY4 plays an important role in cross-presentation in vivo, and is required for antigen processing and presentation and cellular response to virus [40]. Similar to WDFY3, WDFY4 may be involved in the regulation of autophagy, and simultaneously, multiple studies have shown that autophagy was associated with lymphocyte activation, especially in B cells [41]. CD19 is a reliable marker for pre-B cells, and express in most B-cell malignancies [42]. MS4A1/CD20 is a marker gene for B cells and highly expressed in S01 of B cells (Figure 2E), and a recent study had described a subset of CD20+ T cells with immune-regulatory and pro-inflammatory activity [43].

Meanwhile, there are several interleukin-related genes, including IL16, IL2RB, IL18BP, that were significantly highly expressed in S01 of CD8 T cells (Figure 2F). In this cellular state, Macrophage-expressed gene 1 (MPEG1), which was reported to be indispensable for bacterial killing following phagocytosis [44], was also highly expressed (Figure 2F). We next performed functional enrichment analysis on the marker genes of S01 in both B cells and CD8 T cells. We found that they were mainly enriched in the pathways related to immune response (Figure 2G), including response to virus, immune response-activating signal transduction. These results suggested that the marker genes in S01 of both B and CD8 T cells play important roles in MCC.

Fibroblast and monocytes/macrophages cellular states in MCC

Cancer-associated fibroblasts (CAFs) promote cancer development and metastasis [45]. In addition, fibroblasts exhibit a CAF-like phenotype in MCC tumour [46]. Mon/Mac are another important immune cells in TME, which are regarded as tumour-associated macrophages (TAMs) [47]. Emerging studies suggested that TAMs promote tumour growth and subvert anti-tumour immunity [48]. To obtain the phenotypes of fibroblast and Mon/Mac in MCC, we applied EcoTyper to perform in-depth analysis of both cellular types. Two cellular states were identified in fibroblasts (Figure 3A). However, we did not observe an obvious sample concentration trend in cellular states of fibroblasts (Figure 3B). Next, we respectively analyzed the functions of the highly expressed genes in the two cellular states. We found that the highly expressed genes were enriched in similar functions (Figure 3C).

Figure 3. The cellular states landscape of fibroblasts and monocytes/macrophages in MCC. (A) Heatmap depicting fibroblasts cell states identified from MCC bulk transcriptomes. (B) The sankey diagram showed the relationship of fibroblasts cellular states, tissue, and merkel cell polyomavirus status. (C) Barplots showing the enriched functional pathways for genes highly expressed in S01 states and S02 states of fibroblasts. (D) Heatmap depicting Mon/mac cell states identified from MCC bulk transcriptomes. (E) Heatmap showing the expression of the top the most prevalent genes in three Mon/mac cell states.

Moreover, we identified three cellular states of Mon/Mac in MCC (Figure 3D). We found different genes that were highly expressed in each cellular state. Thus, we investigated the expression of top ten ranked genes in cellular states of Mon/Mac in detail. We found that the marker genes showed higher expression in the corresponding cellular states (Figure 3E). In addition, we found that SLIT2, NR1H3, PLA2G7, FGL2 were significantly highly expressed in S02 (Figure 3F). SLIT2 plays an antitumor role in a variety of tumours [49] and recent study has shown that SLIT2 regulates antitumor immunity by altering bone marrow-derived macrophages (BMDMs) [50]. PLA2G7, an enzyme with anti-inflammatory and proinflammatory functions, has been reported to have tumour predictive ability in different tumours [51]. NR1H3 was found to be associated with a subclass of macrophages with better survival [52]. Accumulating evidences have indicated high expression of FGL2 is associated with better prognosis [53]. These results suggested that the subtyping of Mon/Mac can identify cancer-related biomarkers in MCC.

Functional profiling of cellular states in MCC

After identifying cellular states and gene expression analysis for key cellular states, we next tried to predict the potential functions of each cellular state. First, we performed functional enrichment analysis for marker genes of cellular states with hypergeometric test. We found that the majority of genes were enriched in the immune-related pathways and cancer hallmark-related pathways (Figure 4A). In particular, the marker genes of S01 in epithelial cells and PMN cells enriched in cancer-related pathways, including apoptosis, notch signalling, epithelial mesenchymal transition (Figure 4A). Several studies have explored that notch signalling regulates tissue development and homeostasis, and was associated with the occurrence and progression of various malignant tumours [54]. Moreover, different degrees of epithelial-mesenchymal transition (EMT) were observed in Merkel cell lines [55].

Figure 4. Comprehensive functional profiling of cellular states in MCC. (A) Circos plot showing the cancer hallmark and immune-related pathways enriched genes highly expressed in cellular states in MCC. The cell states of each cell type are shown on the left and pathways are on the right, widths of the links indicate the proportion of genes enriched in pathways. (B) Scatter plot showing the normalized enrichment scores (NES) of pathways in epithelial cells S01. Pathways are ranked by NES. Colours of dots are corresponding to the p-values. (C) The enrichment score (ES) distribution for the genes differentially expressed in epithelial cells S01 enriched in five pathways. Each line is for a pathway.

We next conducted an in-depth analysis of epithelial cells S01 and several cancer hallmark-related pathways were found to be significantly enriched in the marker genes of the epithelial S01 (Figure 4B–C), including coagulation, inflammatory response, interferon alpha response. Several oncogenic pathways, including UV response UP, were mainly enriched in genes with low expression in S01 (Figure 4B). As it has been mentioned above that epithelial stem cells may also be the origin cells of MCC, then these results further demonstrated their functions in MCC. Taken together, the comprehensive functional landscape in MCC was revealed by functional enrichment analysis of marker genes in each cellular state.

Discovery of ecosystems in MCC

TME is a complex system composed of heterogeneous cell types, and plays critical roles in cancer development and progression. Therefore, it is essential to uncover the underlying mechanisms of TME in MCC. We employed EcoTyper to decode the multicellular community ecotypes in MCC, which were considered as ecotypes. It turns out that we decoded four ecotypes (Figure 5A). We next analyzed the composition of ecotypes by cell abundance profiles. We found significant differences in the composition of each ecotype, and the majority of MCC patients were divided into specific ecotypes (Figure 5A).

Figure 5. Landscape of cellular ecosystems in MCC. (A) Cellular state abundance patterns of MCC patients, with cell states organized into four ecotypes and tumour samples (columns) ordered by the most abundant ecotype per sample. (B) The proportions of cell types in four ecotypes in MCC. (C) Sankey diagram showing distribution of cellular states in four ecotypes. (D) Network organization of cellular states in four ecotypes. (E) Sankey diagram showing the relationship of four ecotypes, tissue, and merkel cell polyomavirus status. (F) Molecular characteristics of four ecotypes. The abundance of cell types in four ecotypes (top), the functional pathways of four ecotypes (bottom).

To further explore the differences in cell type composition within each ecotype, we calculated the proportion of cell types for per ecotype. We found that both the type and number of cell types in each ecotype were significantly different (Figure 5B–C). For example, E1 consisted of seven cell types including fibroblasts and Mon/Mac, whereas E2 and E3 contained only three cell types (Figure 5C). In addition, the composition of cellular states in ecotypes was also significantly different, with each CE containing 3 to 7 cell states (Figure 5D).

Next, we analyzed the characteristics of patients from different ecotypes, and found that most of the patients in E1 were primary MCC, although there was no obvious difference in viral infection (Figure 5E). To better understand the functions of ecotypes, we performed functional enrichment analysis using hypergeometric tests, including immune- and cancer hallmark-related pathways. We found that genes highly expressed in E1 were enriched in several cancer hallmark-related pathways (Figure 5F, i.e. TNFA signalling via NFKB, inflammatory response). Several genes were also significantly enriched in immune-related pathways, including antimicrobials, cytokines, and cytokine receptors, probably because of the influence of B cells (Figure 5F). At the same time, we found that the high expression genes of E4 were only enriched in the cancer-related pathways. These results suggested that genes in each ecotype play different functional landscapes in MCC.

Transcriptional regulation of ecosystems in MCC

We next explored the function of ecotypes in terms of upstream transcriptional regulators in MCC. To achieve this goal, transcriptional regulation analysis was performed based on marker genes of ecotypes. In total, we found several TF regulators were enriched in two ecotypes (Figure 6(A–B)). E2F transcription factor family members play critical roles in cell cycle regulation and apoptosis [56]. E2F family genes have increased expression levels in a variety of cancers [57]. We found that the target genes of multiple E2F family genes, including E2F1, E2F3 and E2F7, were significantly enriched in E1 highly expressed genes (Figure 6A). SPI1 is a transcriptional activator that may be specifically involved in the differentiation or activation of macrophages or B cells, the targets of SPI1 were highly expressed E3 gene (Figure 6B).

Figure 6. Transcriptional regulators of ecosystems in MCC. (A-B) Representative TFs potentially regulating the genes highly expressed in corresponding ecotypes. Top five TFs were shown for each ecotype of MCC. (A) for E1 and (B) for E3. (C) Heatmap showing the expressions of TFs across patients within different ecotypes. (D) Boxplots showing the expressions of representative TFs. (E-F) Functional enrichment of the target genes of transcription factors. GO terms were clustered into groups based on gene similarities. (E) for E2 and (F) for E3.

Next, we compared the expression of top five TFs ranked by enrichment scores in each ecotype. We found that these TFs were significantly highly expressed in the corresponding ecotypes (Figure 6C). For example, E2F1, E2F3, E2F4, E2F7, TFDP1 have higher expression levels in E2, and SPI1 is significantly highly expressed in E3 (Figure 6D). Moreover, we performed functional enrichment analysis on the target genes of these transcription factors. For E2 ecotype, we found that genes were mainly enriched in DNA replication and related to regulation pathways (Figure 6E). T cell activation and other immune-related pathways were enriched in the target’s genes of SPI1 in E3 (Figure 6F). These results suggested that TFs play important roles in regulating the genes that functions in different ecotypes of TME in MCC.

Discussion

Merkel cell carcinoma (MCC) is a rare aggressive skin cancer associated with immunosuppression [58]. TMEs play key roles in the occurrence, development, invasion and treatment of cancer [59]. The efficacy of cancer immunotherapy largely depends on the TMEs [60]. Comprehensive analysis of MCC based on TMEs can provide new insights into the clinical diagnosis and treatment of MCC. In this study, we stratified MCC patients according to transcriptome-defined subtypes and cell type-specific gene expression patterns in cancer. We used the machine learning method, EcoTyper, to characterize the TMEs in 74 patients of MCC. In total, we comprehensively identified 30 transcriptome-defined cellular states based on 12 cell types, and four multicellular communities. These cellular states and ecosystems provided valuable resources for characterization of TMEs in MCC.

A number of studies have suggested that MCC can be controlled by the immune system, and the markers of different immune cell types may serve as prognostic factors in patients with MCC [61]. We applied EcoTyper to B cells and CD8 T cells, two important immune cell types, and obtained two B cell states and three CD8 T cell states. We found that the majority of patients in S01 of both cell types had better survival of primary MCC and MCPyV + MCC. Moreover, the results of functional enrichment analysis showed that the marker genes of S01 in both B cells and CD8 T cells were significantly enriched in the pathway related to immune response. By in-depth analysis of the two immune cell types, we found that S01 states of both B cells and CD8 T cells play important roles in MCC.

Cancer-associated fibroblasts (CAFs) and tumour-associated macrophages (TAMs) are important components of tumour microenvironment [62,63], affecting cancer progression and metastasis [64,65]. Therefore, EcoTyper was used to analyze fibroblasts and Mon/Mac. We identified two fibroblast cellular states and three Mon/Mac cellular states. We also found that SLIT2, NR1H3, PLA2G7, FGL2 were significantly highly expressed in Mon/Mac cellular states S02. The results suggested that the Mon/Mac subtype may be a potential prognostic factor for MCC. Functional enrichment analysis revealed that genes highly expressed in cellular states were significantly enriched in the cancer hallmark pathways and immune-related pathways.

Moreover, to systematically analyze the TFs of ecotypes, we conducted transcriptional regulation analysis based on marker genes of each ecotype. We obtained some important TFs, including E2F transcription factor family members (e.g. E2F1, E2F3, E2F7) and SPI1. These results suggested that TFs play important roles in the TME of MCC. In this study, we depicted cellular states and ecotypes based on only 74 MCC patients, and more datasets are needed for validation. Moreover, whether the biomarkers identified from B cells, CD8 T cells, fibroblasts and Mon/Mac can be used as prognostic factors needs to be verificated with relevant clinical data.

In conclusion, we conducted comprehensive analysis of TME in MCC based on the transcriptome and a machine learning framework, providing a new perspective for the clinical diagnosis and treatments of MCC.

Authors’ contribution

PZ and JL proposed and managed the study. SC and SL wrote the manuscript. SC and SL collected the material and performed the bioinformatic analysis. RW and SC provided the explanations of the results. Data analysis consisted of SL, SC and PY. All authors contributed to the article and approved the submitted version.

Acknowledgments

We thank the Hainan Province Clinical Medical Center and Academician Workstation in Hainan Province for their support and assistance with this study.

Disclosure statement

No potential conflict of interest was reported by the author(s). The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Data availability statement

Public gene expression profiles used in this work can be acquired from Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/). The datasets used and/or analyzed during the present study are available from the corresponding authors on reasonable request.

Additional information

Funding

This work was supported by National Science Foundation of China (82260474), Hainan Province Science and Technology Special Fund (ZDYF2022SHFZ132).