Multiple datasets to explore the tumor microenvironment of cutaneous squamous cell carcinoma

1.   Introduction

As a non-melanoma skin cancer, cutaneous squamous cell carcinoma(cSCC) has a considerable influence on worldwide health systems. According to linked studies, non-melanoma skin cancer was the most frequent malignancy in men and women in 2017, with up to 7.7 million instances, of which 5.9 million cases were basal cell carcinoma and 1.8 million cases were cSCC ^[1]. In 2020, 1.2 million people develop cSCC, according to a recent study ^[2]. Even though cSCC has a low risk of lymphatic metastasis and a favorable prognosis, 2.1 to 5% of skin cancers advance to a metastatic state ^[3,4]. As the most frequent type of metastatic skin cancer, cSCC has a terrible prognosis at the late stage ^[5], accounting for 20% of cutaneous carcinoma-related deaths ^[4]. CSCC is a type of keratin-forming cell tumor that accounts for around 20% of all keratin-forming cell malignancies ^[6]. Long-term exposure to ultraviolet radiation, irregular circadian rhythms, smoking and alcohol misuse, and being overweight are frequent risk factors for skin cancer ^[7,8]. Additionally, it may develop due to various skin conditions such as scars, burn wounds, discoid erythema, or hypertrophic lichen planus ^[9,10,11,12]. Excision surgery is an excellent therapeutic option for extensive squamous skin cancers ^[13]. Due to its propensity for metastasis and poor prognosis in its advanced stages, we must continue investigating its pathophysiology and developing more effective treatments.

Further research revealed that squamous cell carcinomas, such as those found in the esophagus, the lung, the head and neck, and other locations, are complex tissue composed of various distinct cells ^[14,15,16]. We believe that cSCC is similar in tissue complexity to other types of squamous cell carcinoma. That is, squamous tumor tissue contains cancer cells and a diverse array of other cellular components involved in inflammation, angiogenesis, and the immune response, forming a complex network of cellular regulation ^[17].

The tumor microenvironment is constituted of tumor cells, resident and recruited host cells (cancer-associated stromal cells and immune cells), as well as the secretory products of the corresponding cells (e.g., cytokines and chemokines) and extracellular matrix components (ECM). The tumor microenvironment serves as the "soil" for tumor survival, and tumor growth and metastasis are highly dependent on the surrounding environment, analogous to the "seed-soil" relationship ^[18]. Angiogenesis is required for cancer cells to grow and metastasize in the tumor microenvironment ^[19,20]. The tumor microenvironment contains a diverse array of immune cells, which play a critical role in anticancer immunotherapy. In-depth examination of tumor-infiltrating immune cells can generate novel therapeutic avenues and concepts and produce particular anti-tumor immune responses in specific patient populations ^[21,22]. Thus, tumor cells and the tumor microenvironment complement one another, and the heterogeneity of the tumor microenvironment and the variable response rates of tumor patients to immunotherapy are inextricably linked. We can improve clinical treatment and achieve precision medicine by typing the tumor immune microenvironment.

We now have a better grasp of the molecular mechanisms and expression of oncogenes in the tumor microenvironment because of advances in research methodologies in recent years. Yan's study, for example, used single-cell RNA-sequencing technology(scRNA-seq) to investigate cSCC at the transcriptome level, collecting 350 single-cell transcriptome profiles from skin squamous cancer cells and normal tissue cells to investigate the fundamental biological pathways at play. Members of the S100 gene family, SPRR gene family, and FABP5 gene family are upregulated in cSCC cells, suggesting that S100A9 and FABP5 could be potential therapeutic targets for cSCC ^[23]. Andrew's research is the first to combine geographic mapping of spatial transcriptomics (ST) and multiplexed ion beam imaging (MIBI) with scRNA-seq to create the most comprehensive single-cell transcriptional atlas of cell subpopulations inside cSCC to date, and to analyze the ecological niche in which these cells live. Patients with high expression of tumor-specific keratinocyte (TSK) markers ITGB1 and PLAU had significantly poorer progression-free survival after PD-1 inhibitor treatment, according to this study ^[24].

Given the importance of the immune microenvironment in carcinogenesis and the lack of information about cSCC-associated immune microenvironment, we investigated the interactions between the two and the mechanisms involved. We first investigated the molecular network mechanisms underlying the tumor microenvironment of cSCC using multiple datasets, performing immune cell infiltration analysis, WGCNA analysis, GO and KEGG functional analysis, and developing predictive models for the tumor microenvironment using lasso regression. We used GSEA analysis to compare the route differences between the high and low expression groups. Then we used enrichment analysis and Motif-TF annotation to identify transcription factors that can alter core gene regulation. Additionally, we screened essential genes based on their GO semantic similarity, did a differential analysis on immune-regulated genes, and investigated the relationship between core genes and immune regulation.

2.   Material and methods

2.1.   Data download

We retrieved the GSE42677 Series Matrix File data file from the NCBI GEO public database, along with the annotated file GPL571, and included the expression profile data for 20 patient groups, including the regular group (n = 10) and the tumor group (n = 10). Additionally, we obtained the GSE45164 Series Matrix File data file, annotated with GPL571, which contained the expression profiles of 13 patient groups, comprising normal (n = 3) and tumor (n = 10) groups. To adjust the microarray, we employed the SVA method. Then, we obtained the GSE53462 dataset, which was annotated using the platform GPL10558 and had ten distinct groups of transcriptome data, comprising a regular group (n = 5) and a tumor group (n = 5). We retrieved the GSE45216 dataset for subsequent validation, annotated on GPL570, containing 40 sets of transcriptome data, including AK group (n = 10) and SCC group (n = 30).

2.2.   Immune cell infiltration analysis

The CIBERSORT method is a well-established technique for determining the immune cell types present in the microenvironment. The method is based on support vector regression, and deconvolution is used to assess the expression matrix of immune cell subtypes. It contains 547 biomarkers that distinguish 22 distinct subpopulations of human immune cells, including T-cell, B-cell, plasma cell, and myeloid cell. The CIBERSORT method was used to infer the 22 immune infiltrating cells' relative proportions and perform spearman correlation analysis on gene expression and resistant cell content in this study ^[25].

2.3.   WGCNA analysis

By developing weighted gene co-expression networks, we identified co-expressed gene modules. We investigated the relationship between gene networks and phenotypes, as well as the core genes in the networks. We constructed a co-expression network for all genes in the GSE131761 dataset using the WGCNA-R package and selected genes with the top 5000 variances for further analysis using this technique, with the soft threshold set to 11. To evaluate network connectedness, we turned the weighted adjacency matrix into a topological overlap matrix (TOM) and used hierarchical clustering to generate the TOM matrix of the clustering tree structure. The different branches of the clustering tree correspond to distinct gene modules, and the various colors correspond to specific modules. We identify genes based on their weighted correlation coefficients and put genes with similar expression patterns into a module, splitting tens of thousands of genes into many modules based on their gene expression patterns.

2.4.   GO and KEGG analysis

We used the Metascape database (www.metascape.org) for annotation and visualization, Gene Ontology (GO) analysis, and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis for specific genes to determine the biological activities and signaling pathways involved in illness development. We defined statistical significance as Min overlap ≥ 3.00 or P ≤ 0.01. Both GO/KEGG and Lasso prediction models used genes from the salmon module for 33 samples.

2.5.   Predictive model construction

We picked the WGCNA modular genes most strongly associated with the trait and further developed prediction models using lasso regression. After combining the expression data for each gene, a scoring system was developed for each patient and weighted using the predicted regression coefficients in lasso regression analysis. We utilized the median score value as the cut-off point for the risk score formula and ROC curves to determine the model's predictive accuracy. We created a WGCNA network based on the patients' T-cell characteristics to investigate the critical regulatory genes in SCC. The function "sft$powerEstimate" was used to determine the soft threshold, which was set to 11.

2.6.   GSEA analysis

GSEA analysis ranks genes based on their differential expression between two types of samples. It then determines whether the predetermined genes are enriched at the top or bottom of this ranking table. In this study, we used GSEA to compare the differential KEGG signaling pathway between the high and low expression groups to better understand the molecular mechanisms underlying the core genes in the two groups, with the number of substitutions set to 1000 and the type of substitution set to phenotype. The software used for GSEA analysis is GSEA4.0, and the background gene set is KEGG v7.5.symbol.gmt. GSEA analysis used all genes from the expression profile.

2.7.   Regulatory network analysis of essential genes

In eukaryotes, transcription initiation is a highly complex process that frequently involves the cooperation of numerous protein components. Together, transcription factors and RNA polymerase II create a transcription initiation complex involved in the transcription initiation process. Transcription factors are classified into two types based on their mode of action. The first category includes universal transcription factors, which form a transcription initiation complex with RNA polymerase II before transcription can begin correctly. Cis-acting elements are sequences found in the flanking regions of genes that can influence gene expression. Cis-acting elements, which regulate gene expression, include promoters, enhancers, regulatory sequences, and inducible regions. The cis-acting element does not encode any protein. It merely supplies an action site for the trans-acting factor to bind with. The current study used the R package cisTarget, with the following parameters, rcistarget.hg19. motifdb. cisbpont.500bp, for the Gene-motif ranking database.

2.8.   GO semantic similarity

We classified proteins by their functional similarity to their interacting proteins based on the similarity of GO semantics utilized for gene annotation. Correlation with gene expression profiles validated GO semantic similarity ^[26], establishing a basis for functional comparison of gene products. As a result, GO semantic similarity is widely used in bioinformatics applications such as protein-protein interaction analysis ^[27], pathway analysis ^[28], and gene function prediction ^[29]. In this section, we quantify the functional similarity of proteins. The term "functional similarity" refers to the geometric mean of the semantic similarity of GO terms in terms of cellular components (CC) and biological pathways (BP). It quantifies the strength of the relationship between each protein and its interacting proteins by considering both function and path. The GOSemSim software was used to determine the semantic similarity of interacting histones in CC and BP ^[30], which assumes the GO topology to perform more accurately ^[31]. The geometric mean of semantic similarity in CC and BP was used to evaluate functional similarity further.

2.9.   Immunohistochemistry of key genes

To further validate our findings, we obtained immunohistochemical staining results of four key genes DTYMK, CDCA8, PTTG1 and MAD2L1 in normal skin tissue and tumor tissue on the Human Protein Atlas (https://www.proteinatlas.org/) website, enabling our core genes to be validated again after model validation.

2.10.   Statistical analysis

The R programming language was used to conduct all statistical analyses in this study (version 4.0). All statistical tests were two-sided, and a significance level of P < 0.05 was considered statistically significant. The present analysis specifically involved two groups of patients (control and tumor groups) and did not perform multiple comparisons, therefore no corrected P values were required.

3.   Results

We retrieved the GEO database's SCC-related datasets GSE42677 and GSE45164 and included 33 patients, 13 in the regular group and 20 in the malignant group. We use the SVA algorithm to fix the chips. The results indicated that the batch effect between chips was abolished when the SVA algorithm corrected for it (Figure 1A, B). The microenvironment mainly consists of fibroblasts, immune cells, extracellular matrix, numerous growth factors, inflammatory factors, and unique physicochemical properties, among others. The microenvironment significantly impacts the disease's diagnosis, survival outcome, and clinical treatment sensitivity. The association between the SCC dataset and immune infiltration was studied to understand better the underlying mechanisms affecting tumor progression. We depict a bar graph of the various types of immune cells found in each sample (Figure 2A). We created a heat map of immune cell correlations, highlighting positive correlations with red and negative correlations with blue, with darker colors suggesting stronger correlations (Figure 2B). Additionally, we discovered that when we compared the significance of microenvironment scores between tumor and regular groups, we found that B cells memory (P = 0.016), T cells regulatory (Tregs) (P = 0.001), Mast cells resting (P = 0.002), Mast cells activated (P = 0.02), and several other immune microenvironment factors were significantly different between groups (Figure 2C).

Figure 1.  Removing batch effects between chips. (A) pre-corrected box line plots of GSE42677 and GSE45164; (B) corrected box line plots of GSE42677 and GSE45164. Blue group representatives GSE42677, red group representatives GSE45164.

DownLoad: Full-Size Img PowerPoint

Figure 2.  Microenvironmental immune infiltration analysis. (A) bar graph of different types of immune cells in each sample; (B) heat map of immune cell correlation; (C) differential analysis of infiltration of different kinds of immune cells in the tumor and regular groups.

DownLoad: Full-Size Img PowerPoint

We next created a WGCNA network and identified gene modules using the tom matrix, resulting in the identification of nine gene modules, including black (253), blue (482), greenyellow (865), grey (304), magenta (208) and purple (865). Modules are available in greenyellow (865), grey (304), magenta (208), purple (166), red (2502), salmon (102) and tan (118). We discovered the strongest association between the salmon module and SCC phenotypes by module-trait analysis (cor = 0.54, P = 0.001) (Figure 3). We performed enrichment analysis using the metascape database. We discovered that genes mainly involved in the following pathways, including mitotic cell cycle process、microtubule cytoskeleton organization involved in mitosis、centrosome、Cell cycle、positive regulation of cell cycle (the top half of Figure 4A). The bottom half of Figure 4A is a cluster network of enriched pathways, where nodes sharing the same cluster are usually close to each other. The interplay of the salmon module's genes is seen in Figure 4B. The interaction network in Figure 4B was obtained from the string database to obtain the interactions between proteins and then visualized by cytoscape.

Figure 3.  (A) A co-expression matrix was constructed to group genes into different modules based on their adjacency and similarity, and this was used to obtain a systematic clustering tree among genes; (B) Based on the tom matrix to detect gene modules, a total of nine gene modules were detected, namely black, blue, greenyellow, grey, magenta, purple, red, salmon, and tan. further exploring the correlation analysis between modules and traits, it was found that the salmon module had the highest correlation with the SCC phenotype (cor = 0.54, P = 0.001).

DownLoad: Full-Size Img PowerPoint

Figure 4.  (A) the top half of Figure 4A: enrichment analysis by metascape database showed that genes were mainly enriched in mitotic cell cycle process, microtubule cytoskeleton organization involved in mitosis, centrosome, cell cycle, positive regulation of cell cycle, cell cycle G2/M phase transition, mitotic metaphase plate congression, DNA conformation change, meiotic nuclear division and other pathways; the bottom half of Figure 4A: cluster network of enriched pathways, where nodes sharing the same cluster are usually close to each other; (B) the interaction network which was obtained from the string database to obtain the interactions between proteins and then visualized by cytoscape.

DownLoad: Full-Size Img PowerPoint

To further characterize the salmon module gene set's essential genes, we employed the datasets GSE42677 and GSE45164 as training sets and GSE53462 and GSE45216 as validation sets. We used lasso regression to do feature screening. The results indicated that lasso regression revealed four genes as hallmark genes of SCC for further study. These four genes were DTYMK, CDCA8, PTTG1, and MAD2L1 (Figure 5A–C). Our study used the lasso algorithm to construct the prediction model. The results indicated that the model built using the four genes had a high diagnostic efficacy, with an area under the AUC curve of 0.9846. We subsequently tested the diagnostic model using the validation set, and the findings indicated that the model was relatively stable, with an area under the AUC curve of 0.8567 and 1, respectively (Figure 6A–C). Since the area under AUC curve greater than 0.8 can also be considered to have a relatively strong predictive power, our results suggest that the prediction model can predict not only the difference between normal vs SCC, but also AK vs SCC.

Figure 5.  Screening of SCC core genes. (A) Ten-fold cross-validation of tuning parameter selection in the LASSO model. (B) LASSO coefficient distribution of differential genes. (C) Coefficients of the four Lasso (core) genes.

DownLoad: Full-Size Img PowerPoint

Figure 6.  (A) constructing ROC curves of the prediction model using the training set (GSE42677 and GSE45164), area under the AUC curve is 0.9846; (B) further validating the ROC curves of the diagnosis model using the validation set GSE45216, area under the AUC curve is 0.8567; (C) further validating the ROC curves of the diagnosis model using the validation set GSE53462, area under the AUC curve is 1.

DownLoad: Full-Size Img PowerPoint

We then examined the individual signaling pathways enriched by the four core genes to determine the molecular mechanisms affecting SCC progression. GSEA analysis revealed considerable enrichment in some relevant pathways. The pathways with the highest CDCA8 expression levels were FRUCTOSE AND MANNOSE METABOLISM and CITRATE CYCLE TCA CYCLE. The paths with the highest levels of DTYMK expression were HUNTINGTONS DISEASE and KEGG PYRIMIDINE METABOLISM. MAD2L1 enrichment pathways include RNA DEGRADATION and CELL CYCLE. PTTG1 enrichment pathways are highly expressed in PARKINSONS DISEAS and PATHOGENIC ESCHERICHIA COLI INFECTION (Figure 7).

Figure 7.  Gene set enrichment analysis of four core genes. (A) highly expressed CDCA8-enriched pathway (KEGG FRUCTOSE AND MANNOSE METABOLISM) (B) highly expressed CDCA8-enriched pathway (KEGG CITRATE CYCLE TCA CYCLE); (C) highly expressed DTYMK-enriched pathway (KEGG HUNTINGTONS DISEASE); (D) pathway highly enriched in DTYMK (KEGG PYRIMIDINE METABOLISM); (E) pathway highly enriched in MAD2L1 (KEGG RNA DEGRADATION); (F) pathway highly enriched in MAD2L1 (KEGG CELL CYCLE); (G) pathway enriched by high expression of PTTG1 (KEGG PARKINSONS DISEAS); (H) pathway enriched by high expression of PTTG1 (KEGG PATHOGENIC ESCHERICHIA COLI INFECTION).

DownLoad: Full-Size Img PowerPoint

We analyzed four essential genes and discovered that they are regulated by various transcription factors and other standard processes. As a result, we used cumulative recovery curves, Motif-TF annotation, and a selection of essential genes to do enrichment analysis on these transcription factors. The research revealed that RORA, RORB, and RORC were the primary regulators of the gene set, with MOTIF being annotated as cisbp M5032. We identified two model genes with a normalized enrichment score (NES) of 7.9 enriched in this motif. We display all of the motifs and transcription factors associated with the model genes (Figure 8). Additionally, we screened essential genes based on their GO semantic similarity, and the results indicated that CDCA8 and PTTG1 were more critical across the network (Figure 9).

Figure 8.  Analysis of Motif transcriptional regulation of core genes. (A, B) the red line is the mean of the recovery curve of each motif, the green line is the mean ± standard deviation (sd), and the blue line is the recovery curve of the current motif. The maximum distance point (mean ± sd) between the current motif and the green curve is the maximum enrichment level selected; (C) four core genes enriched to motifs. NES: normalized enrichment score of motifs in the gene set, AUC: area under the curve (used to calculate NES).

DownLoad: Full-Size Img PowerPoint

Figure 9.  Screening of critical genes by GO semantic similarity showed that CDCA8 and PTTG1 are critical in the whole network. The upper and lower limits of the boxes show the 75th and 25th percentiles, and the line in the middle of the box indicates the mean value of similarity. The top two proteins are considered to be key proteins.

DownLoad: Full-Size Img PowerPoint

We did a differential study of immune regulatory genes and discovered that HLA-DMA, HLA-DOA, HLA-DOB, and HLA-DQB2 expression levels were significantly different between the two patient groups (Figure 10A). To investigate the relationship between essential genes and immune regulation, we used Pearson correlation analysis. The results indicated that PTTG1 and HLA-DMA were considerably favorably correlated (Pearson r = 0.62), whereas CDCA8 was strongly negatively connected with HLA-DQB2 (Pearson r = −0.63) (Figure 10B). We used the Human Protein Atlas website (https://www.proteinatlas.org/) to download immunohistochemical staining tissue images of the four core genes in normal skin tissue and tumor tissue, from which we can see the significant differential expression of the four core genes between normal and tumor tissues, further validating the accuracy of our model (Figure 11A–D).

Figure 10.  Correlation analysis of core genes and immune-related genes. (A) Differential analysis of immune regulatory genes showed that HLA-DMA, HLA-DOA, HLA-DOB and HLA-DQB2 were significantly different in the two groups of patients. * denotes P < 0.05, ** denotes P < 0.01, *** denotes P < 0.001; (B) Correlation analysis of critical genes PTTG1 and CDCA8 showed that PTTG1 and HLA-DMA were significantly positively correlated, and CDCA8 was significantly negatively correlated with HLA-DQB2.

DownLoad: Full-Size Img PowerPoint

Figure 11.  Immunohistochemical expression results of 4 core genes from Human Protein Atlas (https://www.proteinatlas.org/). (A) CDCA8 expression in the normal and tumor tissue (HPA028058); (B) DTYMK expression in the normal and tumor tissue (HPA042593); (C) MAD2L1 expression in in the normal and tumor tissue (HPA003348); (D) PTTG1 expression in the normal and tumor tissue (HPA045034).

DownLoad: Full-Size Img PowerPoint

4.   Discussion

Epithelial tumors account for 90% of human malignancies, and cSCC serves as the prototype of epithelial tumors due to its basic features of tissue polarity destruction and basement membrane infiltration. CSCC, the second most common cancer in the United States, affects more than one million people each year ^[32]. CSCC has a decent prognosis in its early stages, but about 4% of patients acquire lymph node metastases, and 1.5 percent of patients with advanced metastases die ^[3]. Despite the efficacy of chemotherapy, immunotherapy, and surgery in treating cSCC, its prognosis remains dismal at the moment ^[33]. Thus, it is critical to understand the tumor microenvironment and pathophysiology of cSCC to enhance the disease's prognosis and attain precision medicine.

We used the GEO dataset on cutaneous squamous cell cancer in this investigation. The study included 13 patients from the control group and 20 patients from the tumor group after batch effects between datasets were reduced using the Combat algorithm. We used the CIBERSORT method to measure immune cells in cSCC and built a WGCNA network based on patient T-cell correlations, and we discovered that the salmon module was the most significantly associated with T cells after evaluating nine gene modules and features. We developed a prediction model based on these four core genes shortly after using lasso regression to screen the modules of salmon for traits and discovered that DTYMK, CDCA8, PTTG1 and MAD2L1 play key roles in cSCC development and progression. The ROC curves of the training and validation sets showed good model stability. We studied the specific signaling pathways enriched in the four core genes using GSEA analysis to infer illness progression biological processes. RORA, RORB and RORC were identified as significant transcription factors after enrichment analysis and Motif-TF annotation of transcription factors regulating the core genes. Following that, we ran a GO semantic similarity screen on them and discovered that CDCA8 and PTTG1 play important roles. We also looked at the association between core genes and immunomodulation and discovered that PTTG1 was positively correlated with HLA-DMA, whereas CDCA8 was dramatically negatively correlated with HLA-DQB2. Finally, we double-checked our findings using the HPA (Human Protein Atlas) online database.

We analyzed immune cell infiltration in the tumor microenvironment of cSCC. We discovered that there was a significant difference between the tumor and regular groups in B cells memory (P = 0.016), T cells regulatory (Tregs) (P = 0.001), Mast cells resting (P = 0.002), Mast cells activated (P = 0.02), and several other immune microenvironment factors. All of these immune cells play a relatively critical role in the development of cSCC. Most cancers, including cSCC, contain Treg T cells, and their presence usually indicates a bad prognosis ^[34,35]. Treg T cells can inhibit cytokine release and proliferation in effector cells. They can also suppress the activation of effectors by influencing antigen-presenting cells. Skin infiltration can be induced by UV irradiation, which causes infiltrating T cells to generate Treg T cells. TGF-β is produced by malignant cells during the carcinogenesis of cSCC and can cause infiltrating effector T cells to become Treg T cells ^[36,37]. Treg T cells were discovered in greater numbers in cSCC than in normal skin, and in moderately as well as poorly differentiated cSCC than in highly differentiated cSCC ^[36]. Macrophages are the most common type of leukocyte found in most cancers. Tumor-associated macrophages (TAMs) are altered by paracrine signals such as TGF-β and IFN-γ in late stages of cSCC carcinogenesis, resulting in inefficient phagocytosis and antigen presentation, and hence acquire an immunoregulatory phenotype ^[38,39]. Around and within the cut margins of surgically excised cSCCs, we discovered a substantial number of M1 and M2 macrophage phenotypes ^[40,41]. TAMs directly stimulate tumor growth by secreting pro-vascular growth factors like vascular endothelial growth factor (VEGF) and matrix metalloproteases (MMP), which breach the normal tissue barrier, allowing for additional invasion ^[40,42]. Because of the ability to create antibodies and effector cytokines, mediate antigen presentation, and contribute to T cell polarization, B cells perform an anti-tumor role ^[43,44]. Furthermore, B cells are involved in tumor growth and cancer immunity ^[45]. B cell depletion causes oncogenic failure in animal models of cSCC, owing to their synthesis of immunoglobulins with chemotactic and pro-inflammatory effects, as well as the deposition of stromal immune complexes early in tumor growth, leading to myeloid activation via Fc receptors ^[46,47]. TNF-α and IL-10 can also be used by B cells to inhibit the immunological response to cSCC ^[48]. Furthermore, it has been discovered that during chemical carcinogenesis, B cells can block the development of cSCC by generating IgE ^[49].

We discovered that DTYMK, CDCA8, PTTG1 and MAD2L1 play critical roles in cSCC development and progression. DTYMK (deoxythymidylate kinas) may be used to predict the prognosis of non-small cell lung cancer and hepatocellular carcinoma patients. It is closely associated with cell cycle and acid metabolism pathways and suppresses the immune microenvironment during tumorigenesis. Additionally, overexpression of DTYMK in some patients increases their sensitivity to chemotherapeutic agents, suggesting that it may be a potential therapeutic target ^[50,51]. CDCA8 (cell division cycle associated 8) may provide a novel molecular method for the initial treatment of hepatocellular carcinoma and preventing metastases. CDCA8 silencing reduces the growth of hepatocellular carcinoma via restoring the ATF3 oncogene and inactivating oncogenic AKT/β-linked protein signaling, whereas CDCA8 overexpression promotes tumor cell proliferation. Additionally, CDCA8 is a critical regulator of mitosis and is being investigated as a possible predictive biomarker for many malignancies, including lung, breast, and colon ^[52]. PTTG1 (pituitary tumor-transforming gene-1) is highly expressed in esophageal, laryngeal, lung, and hepatocellular carcinomas ^[53,54], and is associated with the degree of differentiation, invasion, and metastasis of the tumor ^[55,56]. PTTG1 is a chromosomal stability regulator that can aggravate tumor development progression, and inhibiting PTTG1 decreases cancer cell proliferation and promotes apoptosis. Additionally, PTTG1 expression is increased in hematologic malignant proliferative disorders such as multiple myeloma and is significantly related to a bad prognosis ^[57]. MAD2L1 (mitotic arrest deficient 2 like 1) is an emerging target for anticancer therapy because it plays a critical role in mitosis, the cell cycle checkpoint, and the DNA damage response ^[58]. MAD2L1 has been demonstrated to be a novel prognostic factor for lung adenocarcinoma and squamous cell carcinoma, and it also interacts with CDC20 and BUB1B ^[59,60,61]. MAD2L1 and MAD2L2 can result in aberrant salivary ductal cancer formation ^[62,63].

In summary, our current work elucidates the immune cell types and interactions in the tumor microenvironment of cSCC, identifies critical components, and identifies enrichment routes via core genes. Meanwhile, we conducted differential analysis on immune regulatory genes and connected them with disease-associated genes.

5.   Conclusions

In summary, we discovered that DTYMK, CDCA8, PTTG1 and MAD2L1 are the critical genes in cSCC, with CDCA8 and PTTG1 being particularly important. The primary regulators are RORA, RORB, and RORC. B cells memory, T cells regulatory, Mast cells resting, and Mast cells activation were significantly different between the tumor and regular groups. We compared essential genes with immune regulatory genes and discovered a strong positive correlation between PTTG1 and HLA-DMA and a significant negative correlation between CDCA8 and HLA-DQB2.

Conflicts of interest

All authors declare no conflicts of interest in this paper.