Integrative Analysis of Biomarkers and Mechanisms in Adamantinomatous Craniopharyngioma

Craniopharyngioma is a benign tumor, and the predominant treatment methods are surgical resection and radiotherapy. However, both treatments may lead to complex complications, seriously affecting patients’ survival rate and quality of life. Adamantinomatous craniopharyngioma (ACP), as one of the histological subtypes of craniopharyngioma, is associated with a high incidence and poor prognosis, and there is a gap in the targeted therapy of immune-related genes for ACP. In this study, two gene expression profiles of ACP, namely GSE68015 and GSE94349, were downloaded from the Gene Expression Omnibus (GEO) database. The differentially expressed genes (DEGs) were identified by the Limma package, and 271 differentially expressed immune-related genes (DEIRGs) were obtained from the Immport database. The gene ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and Gene Set Enrichment Analysis (GSEA) were performed for annotation, visualization, and integrated discovery. Five hub genes, including CXCL6, CXCL10, CXCL11, CXCL13, and SAA1, were screened out through protein-protein interaction (PPI) network interaction construction. Two diagnostic markers, namely S100A2 and SDC1 (both of which have the Area Under Curve value of 1), were screened by the machine learning algorithm. CIBERSORT analysis showed that M2 macrophages, activated NK cells, and gamma delta T cells had higher abundance in ACP infiltration, while CD8+ T cells, regulatory T cells, and Neutrophils had less abundance in ACP infiltration. The expression of gamma delta T cells was positively correlated with CXCL6, S100A2, SDC1, and SAA1, while CD8+ T cells expression was negatively correlated with CXCL6, S100A2, SDC1, and CXCL10. ACP with high CXCL6 showed remarkable drug sensitivity to Pentostatin and Wortmannin via CellMiner database analysis. Our results deepened the understanding of the molecular immune mechanism in ACP and provided potential biomarkers for the precisely targeted therapy for ACP.

Introduction

Craniopharyngioma (CP) is a rare epithelial neoplasm derived from the remnants of Rathke’s pouch (Prabhu and Brown, 2005). Due to its proximity to critical neurovascular structures such as the optic nerve, pituitary, hypothalamus, and internal carotid artery, pre-and postoperative complications such as visual impairment, endocrine deficiencies, and hypothalamic dysfunctions frequently occur after CP, thus leading to a low survival rate and poor quality of life in CP patients (Müller et al., 2019).

Currently, the treatment of CP mainly includes gross total resection (GTR) and subtotal resection (STR) followed by adjuvant radiotherapy (RT), but it is controversial which of these two methods is more advantageous. Although GTR reduces the possibility of recurrence to the greatest extent, it is associated with a higher risk of neurological complications such as hypothalamus injury. By contrast, STR + RT may retain patients’ neurological function but is associated with a higher risk of recurrence and radiotherapy complications (Aggarwal et al., 2013; Shi X.e. et al., 2017). Therefore, a novel and rational treatment safer than the current treatment methods is needed for this challenging tumor.

The histological subtypes of CP include papillary craniopharyngioma (PCP) and adamantinomatous craniopharyngioma (ACP), and the incidence of ACP is higher than that of PCP. More than 90% of PCP carries BRAF mutations and BRAF V600E protein expression (Malgulwar et al., 2017). Molecular-targeted therapy with BRAF and MEK inhibitors for PCP with BRAF mutation may achieve the goal of reducing tumor volume, thus becoming a neoteric treatment for PCP (Brastianos et al., 2014, 2016; Roque and Odia, 2017; Juratli et al., 2019). By contrast, studies have reported that activation of the Wnt pathway by CTNNB1 mutations and positive β-catenin protein expression is the pathogenic mechanism of ACP, but there is still no clear molecular therapeutic target for ACP at present.

Immune-related genes play crucial roles in the tumor immune microenvironment (Li et al., 2020), and immunotherapy is a validated and critically important approach for treating patients with cancer (Hegde and Chen, 2020). Coy et al. found that PD-1/PD-L1 played an important role in the pathogenesis of ACP and suggested that PD-1/PD-L1 could be a potential immunotherapy target (Coy et al., 2018); however, the drug sensitivity of this target was not tested in their study. There are few studies on ACP immune-related differential genes, so it is necessary to study further the expression and mechanism of ACP immune-related differential genes to find feasible immunotherapy targets.

In this study, we downloaded two ACP gene expression profile datasets from Gene Expression Omnibus (GEO) database and screened out the differentially expressed immune-related genes (DEIRGs). We analyzed the biological functions of DEIRGs by enrichment analysis and used protein-protein interaction (PPI) network interaction and machine learning algorithm to screen hub genes and diagnostic markers. We calculated the abundance of immune infiltrating cells in tumor cells and analyzed the correlation between the immune infiltrating cells and the biomarkers. In addition, we used the CellMiner database to predict the drug sensitivity of the screened biomarkers. Our findings on the therapeutic targets based on the molecular immune mechanism could be further employed as promising biomarkers or treatment targets for ACP.

Materials and Methods

Data Collection and Preprocessing

ACP expression profiles datasets GSE68015 (Gump et al., 2015) and GSE94349 were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) (Barrett et al., 2013). The species selected were Homo sapiens, and the platform was GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array. The GSE68015 dataset included 15 ACP samples and 16 control brain tissues, and the GSE94349 dataset included 9 ACP samples and 17 control brain tissues. The original data of GSE68015 and GSE94349 datasets were read through the DealGPL570 package, background correction and data normalization were carried out, and the gene expression matrices of the two datasets were obtained, respectively.

Identification of DEIRGs

Differentially expressed genes (DEGs) in the GSE68015 and GSE94349 datasets were screened by the “Limma” package (Ritchie et al., 2015), and volcano maps of DEGs were performed to explore the differentially expressed DEGs using the “ggplot2” package. DEGs satisfy p-value < 0.05 and |log2FC| >1(The cut-off criteria for statistical significance were a log fold change (FC) of greater than one and a p value of less than 0.05). The immune-related genes (IRGs) list was downloaded from the ImmPort database (https://immport.niaid.nih.gov) (Bhattacharya et al., 2018), and then the DEIRGs were extracted from DEGs.

Functional and Pathway Enrichment Analysis

We used the “clusterProfiler” package (Yu et al., 2012) to perform Gene Ontology (GO) (Ashburner et al., 2000) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al., 2017) enrichment analysis on DEIRGs. GO analysis included annotation of biological processes (BPs), Molecular functions (MFs), and Cellular Components (CCs). Also, Gene Set Enrichment Analysis (GSEA) was performed on the gene expression matrix through the “clusterProfiler” package (Hänzelmann et al., 2013), and “c2. cp.kegg.v7.0. symbols.gmt” was selected as a reference gene set. A false discovery rate (FDR) < 0.25 and p < 0.05 was considered to be significantly enriched.

PPI Network Construction and Identification of Hub Genes

These DEIRGs were put into the STRING database (https://www.string-db.org/) (von Mering et al., 2003) to obtain the interaction relationships, and Cytoscape software (Shannon et al., 2003) was used to visualize the interaction relationships of DEIRGs with interaction scores greater than or equal to 0.9. The hub genes in the network were mined based on Molecular Complex Detection (Bandettini et al., 2012) and cytoHubba (Chin et al., 2014).

Screening and Validation of Diagnostic Markers

We used two machine learning algorithms, namely Random Forest and least absolute shrinkage and selection operator (LASSO) regression, to screen diagnostic markers in the GSE68015 dataset. The “glmnet” package (Engebretsen and Bohlin, 2019) was used to implement the LASSO algorithm, and the “randomForest” package (Alderden et al., 2018) was used to implement the Random Forest algorithm. Finally, the intersection of diagnostic markers obtained by the two algorithms was taken as the final result, and receiver operating characteristic (ROC) analysis was performed for each diagnostic marker in the GSE94349 dataset to verify its accuracy.

Relationship Between Immune Cell Infiltration and Diagnostic Markers

We uploaded the gene expression matrix data to CIBERSORT to obtain the immune cell infiltration matrix. A heatmap was made using the R language “pheatmap” package (https://CRAN.R-project.org/package=pheatmap) to show the distribution of 22 kinds of immune infiltrating cells in each sample. The “ggplot2” package was used to map the distribution of immune cells in the samples. The “corrplot” package draws a correlation heatmap to visualize the correlation of 22 immune infiltrating cells. The “ggplot2” package was used to create a violin diagram to visualize differences in the infiltration of 22 immune cells. We analyzed the correlation between the selected hub genes, the diagnostic markers, and the immune infiltrating cells. The results were visualized using the “ggplot2” package (Ito and Murphy, 2013).

Drug Sensitivity Analysis of Hub Gene and Diagnostic Markers

Data of genes and drug sensitivity were downloaded from the CellMiner database (https://discover.nci.nih.gov/cellminer/) (Shankavaram et al., 2009), and the top 16 drugs correlated with the hub genes and the diagnostic markers were selected and visualized by “corrplot” package.

Results

Data Preprocessing and DEIRGs Screening

The overview of the workflow is shown in Figure 1. We extracted 3,727 DEGs from the GSE68015, including 1,421 up-regulated genes and 2,306 down-regulated genes, and 4,242 DEGs from the GSE94349, including 1848 up-regulated genes and 2,394 down-regulated genes. The results of the DEGs analysis are shown in Figure 2. The intersection of DEGs and immune-related genes in these two datasets was taken to obtain 271 DEIRGs.

FIGURE 1

FIGURE 1. Flow chart of methodologies applied in this study.

FIGURE 2

FIGURE 2. Analysis results of DEGs in ACP. Heatmaps of DEGs identified in GSE68015 (A) and GSE94349 (B). Volcano plot analysis identifies DEGs in GSE68015 (C) and GSE94349 (D).

Functional Correlation Analysis

The significance of each gene was evaluated by biological process (BP), molecular function (MF), and cellular component (CC) semantics based on GO analysis. GO analysis showed that DEIRGs were mainly associated with “leukocyte migration” and “cell chemotaxis” (Figure 3A). KEGG analysis results are shown in Figure 3F. The pathways enriched by DEIRGs mainly include “cytokine receptor interaction” and “TNF signaling pathway.” Detailed information of GO and KEGG enrichment analysis is shown in Table.1. GSEA enrichment pathways are mainly involved in the Blalock Alzheimer’s Disease Up and Chen Metabolic Syndrom Network pathways (Figure 4), and GSEA details are shown in Table.2.

FIGURE 3

FIGURE 3. Functional correlation analysis of DEIRGs. Enrichment analysis of GO biological functions (A). The X horizontal axis represents the proportion of DEGs enriched in the GO term. The bluer the color, the greater the corrected p value. The size of the dots represents the number of enriched genes. Other manifestations of GO enrichment analysis (B–E). Enrichment analysis of KEGG pathway (F–H).

TABLE 1

TABLE 1. GO and KEGG pathway enrichment analysis of DEIRGs in the most significant module.

FIGURE 4

FIGURE 4. The results of the GSEA analysis in GSE68015 (A–C) and GSE94349 (D–F).

TABLE 2

TABLE 2. The results of GSEA function enrichment analysis.

PPI Network Construction and Hub Gene Selection

The network was constructed by the interaction of these DEIRGs with nodes with scores greater than or equal to 0.9 in the STRING database. Five hub genes, including CXCL6, CXCL10, CXCL11, CXCL13, and SAA1, were screened out by MCODE and cytoHubba plugins. PPI network interaction and hub gene mining are shown in Figure 5.

FIGURE 5

FIGURE 5. PPI network construction and hub gene mining. (A) PPI network diagram with DEIRGs interaction fraction greater than or equal to 0.9, including 156 nodes and 810 edges. The more purple the color is, the greater the node degree is, and the bluer the color is, the smaller the node degree is. (B) The hub node mined by the MCODE plugin has a total of 23 nodes and 253 edges. (C) hub nodes excavated by cytoHubba plugin on the basis of MCODE. There are altogether 5 nodes and 10 edges. (D) GO enrichment cloud and rain map of 5 hub genes.

We ranked the top 5 genes based on the average functional similarity (Figure 5D), which showed that CXCL11 was the most important gene.

Identification and Validation of Diagnostic Markers

In the GSE60815 dataset, the LASSO regression (Figures 6A,B) and Random Forest (Figures 6C,D) algorithms were used to screen diagnostic markers. Two diagnostic markers were selected by the LASSO method, and 18 were selected by the Random Forest method. The two diagnostic markers, namely S100A2 and SDC1, were selected as the final results (Figure 6E), and the S100A2 and SDC1 genes were analyzed by ROC in the external dataset GSE94349 (Figures 6F,G). The results showed that the Area Under Curve (AUC) values of S100A2 and SDC1 genes were both 1, so we selected S100A2 and SDC1 genes as the final diagnostic markers.

FIGURE 6

FIGURE 6. LASSO and Random Forest algorithms screening for diagnostic markers (A,B). LASSO regression model accuracy curve based on Random Forest (C) and error curve based on Random Forest (D). The Venn diagram shows the intersection of diagnostic markers obtained by the two algorithms (E). The ROC curve of the diagnostic efficacy verification after fitting S100A2 (F) and SDC1 (G) in GSE94349.

Analysis of Immune Cell Infiltration and Its Correlation With Hub Genes Diagnostic Markers

The results of the heat map and bar chart of immune cell infiltration showed significant differences in immune infiltrating cells between ACP and control samples (Figures 7A,C,E,G). Violin diagram on immune cells infiltrating differences (Figures 7D,H) showed that, in ACP samples, M2 macrophages, activated NK cells, and gamma delta T cells infiltration were relatively more. In contrast, there were relatively few CD8+ T cells, regulatory T cells, Neutrophils, and activated mast cells infiltration. The heat map of 22 infiltrating immune cells showed that resting CD4+T cells were negatively correlated with gamma delta T cells and monocytes while positively correlated with resting dendritic cells and M1 macrophages (Figures 7B,F).

FIGURE 7

FIGURE 7. Evaluation and visualization of immune cell infiltration in GSE60815 datasets (A–D) and GSE94349 datasets (E–H). Bar chart of immune infiltrating cells between ACP and control samples (A,E). Correlation heat map of 22 types of immune cells. The size of the colored squares represents the strength of the correlation; blue represents a positive correlation; red represents a negative correlation. The darker the color, the stronger the correlation (B,F). A heat map of immune cells between ACP and control samples (C,G). Violin diagram of the proportion of 22 types of immune cells, and comparison of immune cells between the two ACP and control samples (D,H).

We analyzed and integrated the correlation between the five hub genes and the two diagnostic markers and immune infiltrating cells in two datasets, respectively. The results showed significant differences in various immune cells between the high and low groups with differential hub genes and diagnostic markers (Figures 8, 9). CXCL6 was positively correlated with gamma delta T cells and M0 macrophages but negatively correlated with CD8+ T cells and regulatory T cells; CXCL10 was positively correlated with M1 macrophages but negatively correlated with CD8+ T cells and activated NK cells; CXCL11 was positively correlated with M0 macrophages; CXCL13 was negatively correlated with activated NK cells and Neutrophils; SAA1 was positively correlated with gamma delta T cells and M0 macrophages but negatively correlated with regulatory T cells; S100A2 was positively correlated with gamma delta T cells and M0 macrophages but negatively correlated with CD8+ T cells and regulatory T cells; SDC1 was positively correlated with gamma delta T cells and M0 macrophages but negatively correlated with CD8+ T cells and Neutrophils.

FIGURE 8

FIGURE 8. Correlation analysis of immune infiltrating cells with hub gene and diagnostic markers in GSE60815 dataset. CXCL6 (A), CXCL10 (B), CXCL11 (C), CXCL13 (D), S100A2 (E), SAA1 (F), and SDC1 (G).

FIGURE 9

FIGURE 9. Correlation analysis of immune infiltrating cells with hub gene and diagnostic markers in GSE94349 dataset. CXCL6 (A), CXCL10 (B), CXCL11 (C), CXCL13 (D), S100A2 (E), SAA1 (F), and SDC1 (G).

The Drug Sensitivity Analysis of Hub Genes and Diagnostic Marker

We analyzed the relationship between the five hub genes and the two diagnostic markers and drug sensitivity. It showed that CXCL6 was significantly positively correlated with Pentostatin and Wortmannin, indicating that with the increase of CXCL6 expression (Figure 10), the sensitivity of the human body to Pentostatin and Wortmannin increased.

FIGURE 10

FIGURE 10. CXCL6 expression was positively correlated with Pentostatin (A) and Wortmannin (B) sensitivity.

Discussion

Although CP is a benign tumor, it is associated with a high mortality rate. According to the clinical research of Khafaga Y et al., the mortality rate of CP patients treated with surgery was as high as 16% (Khafaga et al., 1998). Harvey Cushing (Elliott et al., 2011) described CP as “one of the most baffling problems to the neurosurgeon.” CP includes two pathological types: ACP and PCP. The incidence of ACP is significantly higher than that of PCP, and the prognosis of ACP is worse than PCP (Adamson et al., 1990). Conventional surgery and radiation therapy for ACP are associated with complex complications and poor prognosis. This study aims to explore the biological mechanism of ACP and search for meaningful molecular markers, which may help find novel diagnosis methods and safe treatments for ACP.

In this study, we performed an integrative analysis of two expression profile datasets (GSE68015 and GSE94349), including 24 ACP samples and 33 normal samples. 271 DEIRGs in total were identified from these two expression profile datasets. Protein-protein interaction (PPI) network interaction and machine learning algorithm were used to screen out five hub genes and two diagnostic markers. CIBERSORT tools have also facilitated the analysis of immune cell infiltration patterns of ACP, and Gamma delta T cells and CD8+ T cells were found to be closely related to these molecular markers. CXCL6 expression was found positively correlated with the sensitivity of Pentostatin and Wortmannin in the Cellminer database. We sought to explore the potential molecular therapeutic targets and related molecular mechanisms of ACP and find a more optimal therapeutic approach.

A total of 271 DEIRGs were screened out in this study. To investigate the biological processes of the DEIRGs involved in ACP, we conducted enrichment analysis of GO, KEGG, and GSEA, respectively. GO enrichment analysis showed that leukocyte migration and cell chemotaxis signaling pathways were significantly associated with the occurrence and development of ACP. KEGG enrichment analysis showed that cytokine-cytokine receptor interaction and TNF signaling pathway were the significantly enriched pathways. Cytokines include chemokines, tumor necrosis factors (TNF), and interleukins. Chemokines may induce the directed migration of white blood cells, which is involved in inflammation and immune response in tumors (Balkwill, 2004). Gong et al. confirmed that the high expression of various cell chemokines in ACP was associated with poor recurrence-free survival (Gong et al., 2014). TNF-α plays an important role in tumor promotion and progression (Balkwill, 2006). Simone et al. suggested that TNF-α and IL-6 cells and other cytokines synergistically activate STAT3 and NF-κB pathways to promote colorectal cancer cell proliferation (De Simone et al., 2015). Mori et al. confirmed the expression of TNF-α in the cystic fluid of ACP tumors (Mori et al., 2004), which is consistent with our data mining results. Therefore, we believe that cytokines play an essential role in the immune response, inflammation, tumor cell proliferation, and apoptosis in ACP.

We screened out five hub genes, namely CXCL6, CXCL10, CXCL11, CXCL13, and SAA1, by constructing the PPI network. The first four of the five hub genes belong to chemokines, which, as key signaling molecules in the tumor microenvironment, play an important role in the interaction between tumor cells and tumor stromal cells (Balkwill, 2004). CXCL11 and CXCL10 are selective ligands for CXCR3, usually expressed at low levels under homeostasis but up-regulated under cytokine stimulation. CXCL11 and CXCL10 are secreted mainly by monocytes, endothelial cells, fibroblasts, and cancer cells in response to IFN-γ and are synergistically enhanced by IFN-α. Some literature suggests that CXCL11 and CXCL10 can promote the proliferation and metastasis of tumor cells (Cambien et al., 2009; Puchert et al., 2020). It has been reported that tumors with high expression of CXCL10 and CXCR3 simultaneously have stronger invasion and metastasis (Wightman et al., 2015). CXCL10, CXCL11/CXCR3 axis can be used as the target of tumor immunotherapy (Tokunaga et al., 2018). CXCL6, as another chemokine family member, also plays an important role in tumor proliferation and metastasis in various cancer, including non-small cell lung cancer (Li et al., 2018) and colon cancer (Ma et al., 2017). Zheng et al. reported that CXCL6 up-regulates PD-L1 expression by activating the STAT3 pathway, promoting the growth and metastasis of esophageal squamous cell carcinoma cells in vivo and in vitro (Zheng et al., 2021). CXCL13 and its chemokine receptor 5 (CXCR5) are among the key chemotactic factors which play crucial roles in deriving cancer cell biology (Hussain et al., 2019). In addition, the SAA1 gene encodes high-density lipoprotein-associated with inflammation, and high expression of SAA1 in pancreatic cancer predicts poor prognosis [39]. Zhang et al. pointed out that inhibiting the SAA1 gene could down-regulate the AKT signaling pathway and promote glioblastoma cell apoptosis, and SAA1 could be used as a potential therapeutic target and prognostic index of glioblastoma (Tokunaga et al., 2018).

LASSO Logistic and Random Forest Regressions were used to screen molecular diagnostic markers in the GSE68015 dataset, and the GSE94349 dataset was used to verify accuracy. Finally, two molecular diagnostic markers, namely S100A2 and SDC1, were selected. S100A2, a member of the S100 family, is a calcium-binding protein. Studies have confirmed that S100A2 participates in processes such as cell proliferation and migration in some tumors (Bulk et al., 2009; Wolf et al., 2011) A previous study showed that high expression of S100A2 in patients with non-small cell lung cancer has a worse prognosis (Wang et al., 2005). Another study indicated that the high expression of S100A2 is an independent prognostic marker for the poor prognosis of pancreatic cancer (Dreyer et al., 2020). Huang et al. suggested that down-regulating S100A2 was able to inhibit the Wnt/β-Catenin signaling pathway, thus reliving the epithelial-mesenchymal transition (EMT) in pulmonary fibrosis (Huang et al., 2021). The accumulation of β-catenin in the nucleus, which further leads to activation of the Wnt pathway, is the main pathogenic mechanism of ACP (Larkin and Ansorge, 2013). We speculate that S100A2 plays an important role in the pathogenesis of ACP. SDC1 is a type I transmembrane glycoprotein mainly expressed in epithelial cells, and disorder of SDC1 expression leads to cancer development by promoting cell proliferation, metastasis, invasion, and angiogenesis. The expression of SDC1 varies in different types of tumors. SDC1 is low expressed in gastric cancer and colorectal cancer, whereas SDC1 is highly expressed in plasmacytoid urothelial carcinoma, liver cancer, and glioma (Liao et al., 2020). Chen et al. reported that SDC1 overexpression significantly increased glioma cells migration, while SDC1 knockdown had the opposite effects (Chen et al., 2017). Another study confirmed that SDC1 knockdown inhibits glioma cell proliferation and invasion by deregulating a c-src/FAK-associated signaling pathway (Shi S. et al., 2017). Therefore, we believe that SDC1 may play a key role in the pathogenetic process of ACP and may serve as a potential therapeutic target in the treatment of ACP.

In this study, the CIBERSORT is used for the first time to analyze the immune cell infiltration in ACP and normal brain tissue to explore the role of tumor-infiltrating immune cells in ACP. We found that increased infiltration of M2 macrophages, activated NK cells, and a decreased infiltration of regulatory T cells, Neutrophils may be related to the occurrence and development of ACP. Macrophages can be polarised towards either M1macrophages or M2 macrophages, and M2 macrophages play an important role in the proliferation of CNS tumors. Komohara et al. reported that M2 macrophages promote the proliferation of glioma cells through M-CSFR/Stat3 (Komohara et al., 2012). In another study, we found that the expression of M2 macrophages in the recurrent CP samples was higher than that in the primary and recurrent CP samples of the same patient, and the more expression of M2 macrophages in the recurrent CP samples, the shorter the recurrence time of the tumor (Lin et al., 2019). We conclude that M2 macrophages play an important role in ACP tumor proliferation. In this study, the proportion of neutrophils in ACP samples was low and significantly lower than that in the normal control group, which was consistent with the results of Coy S et al. (Coy et al., 2018).

By analyzing the correlation between the hub genes and diagnostic markers and immune cells, it was found that CXCL6, S100A2, SDC1, SAA1 were significantly positively correlated with gamma delta T cells. CXCL6, S100A2, SDC1, CXCL10 were significantly negatively correlated with CD8+ T cells. In recent years, studies have confirmed that gamma delta T cells play an important role in the occurrence and progression of tumors (Fleming et al., 2017), while CD8+ T cells are associated with the prognosis of patients and the outcome of immunotherapy. Daley et al. confirmed that tumor-infiltrating gamma delta T cells promote tumor progression through the PD-1/PD-L1 pathway restraining αβ T cell activation (Daley et al., 2016). Gamma delta T17 cell is the main manufacturer of IL-17. Coffelt et al. confirmed in a breast cancer model that IL-17 might induce tumor immunosuppression and promote tumor neovascularization, among which, Il-17 inhibits the activity of CD8+T cells by regulating the tumor-infiltrating neutrophils (Coffelt et al., 2015). We hypothesized that CXCL6, S100A2, and SDC1 lead to tumor progression and immune escape by raising gamma delta T cells and reducing CD8+T cells, but this hypothesis needs further verification.

Finally, we used CellMiner to analyze the drug sensitivity of hub gene and diagnostic markers, and we found that with the increase of CXCL6 expression, the drug sensitivity to Pentostatin and Wortmannin increased. We believe that in ACP with high CXCL6 expression, Pentostatin and Wortmannin can be used as potential molecular targeted drugs for treatment, but the mechanism of action between drugs and molecules needs to be confirmed by further studies.

With the advent of the era of molecular therapy, the treatment of ACP requires more correlation studies on the precise molecular targeted therapy, especially for complicated and rapidly-recurring ACPs and those not suitable for surgical treatment. Grob et al. found that tocilizumab and bevacizumab could reduce tumor volume in cystic ACP with high IL-6 expression (Grob et al., 2019). Apps et al. reported that inhibition of MAPK/ERK signaling pathway by MEK inhibitors could reduce the proliferation of ACP tumor cells or increase tumor cell apoptosis (Apps et al., 2018). These studies, including this one, provide significant clues and evidence for the feasibility of molecularly targeted therapy for ACP.

There are some limitations in the current study. First, ACP includes the solid component, cystic wall, cystic fluid, calcification, and other pathological tissue components. The expression profile data were not precisely selected from specific tumor components, and the gene expression may differ with different tumor tissue components. Second, this study should apply basic experiments to verify the expression of the excavated hub gene and diagnostic marker gene. For example, the role of hub genes and diagnostic markers should be fully elucidated through real-time PCR, immunohistochemistry, immunofluorescence. Third, due to the lack of clinical characteristics data of ACP patients in the original data mining, it is not possible to carry out the correlation study of hub gene, diagnostic markers, and immune infiltrating cells with clinical characteristics of ACP, such as imaging characteristics, recurrence, and survival rates.

Through integrative analysis of biomarkers and mechanisms in ACP expression profile datasets in this study, we believe that the results may provide new insights into the diagnosis and treatment of ACP. Although some results of the previous studies are consistent with the results of our analysis, further experimental validation on the reliability of the results is needed.

Conclusion

In summary, we found five hub genes, namely CXCL6, CXCL10, CXCL11, CXCL13, and SAA1, and two diagnostic markers, namely S100A2 and SDC1, play a crucial role in the pathogenesis of ACP. We also found that M2 macrophages, activated NK cells, gamma delta T cells, CD8+ T cells, regulatory T cells, Neutrophils, and activated mast cells were involved in the occurrence and progression of ACP. Among them, gamma delta T cells were significantly positively correlated with CXCL6, S100A2, SDC1, and SAA1, while CD8+ T cells were significantly negatively correlated with CXCL6, S100A2, SDC1, and CXCL10. Moreover, ACP with high CXCL6 expression showed remarkable drug sensitivity to Pentostatin and Wortmannin; thus, we believe that Pentostatin and Wortmannin can be used as potential molecular targeted drugs to treat ACP with high CXCL6 expression, but the specific drug mechanism needs further verification. The hub genes, diagnostic markers, and immune cells identified in this study may serve as targets for molecular therapy or immunotherapy and provide guidance for developing a novel and safe treatment for ACP.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Author Contributions

HZ, DL, and WZ designed the study and interpreted data. DL, WZ, JY, and HW conducted bioinformatic and statistical analyses. DL and WZ wrote the paper. All authors reviewed and approved the submitted version.

Conflict of Interest

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.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Adamson, T. E., Wiestler, O. D., Kleihues, P., and Yaşargil, M. G. (1990). Correlation of Clinical and Pathological Features in Surgically Treated Craniopharyngiomas. J. Neurosurg. 73, 12–17. doi:10.3171/jns.1990.73.1.0012

PubMed Abstract | CrossRef Full Text | Google Scholar

Aggarwal, A., Fersht, N., and Brada, M. (2013). Radiotherapy for Craniopharyngioma. Pituitary 16, 26–33. doi:10.1007/s11102-012-0429-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Alderden, J., Pepper, G. A., Wilson, A., Whitney, J. D., Richardson, S., Butcher, R., et al. (2018). Predicting Pressure Injury in Critical Care Patients: A Machine-Learning Model. Am. J. Crit. Care 27, 461–468. doi:10.4037/ajcc2018525

PubMed Abstract | CrossRef Full Text | Google Scholar

Apps, J. R., Carreno, G., Gonzalez-Meljem, J. M., Haston, S., Guiho, R., Cooper, J. E., et al. (2018). Tumour Compartment Transcriptomics Demonstrates the Activation of Inflammatory and Odontogenic Programmes in Human Adamantinomatous Craniopharyngioma and Identifies the MAPK/ERK Pathway as a Novel Therapeutic Target. Acta Neuropathol. 135, 757–777. doi:10.1007/s00401-018-1830-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Ashburner, M., Ball, C. A., Blake, J. A., Botstein, D., Butler, H., Cherry, J. M., et al. (2000). Gene Ontology: Tool for the Unification of Biology. Nat. Genet. 25, 25–29. doi:10.1038/75556.Gene

PubMed Abstract | CrossRef Full Text | Google Scholar

Balkwill, F. (2004). Cancer and the Chemokine Network. Nat. Rev. Cancer 4, 540–550. doi:10.1038/nrc1388

PubMed Abstract | CrossRef Full Text | Google Scholar

Balkwill, F. (2006). TNF-α in Promotion and Progression of Cancer. Cancer Metastasis Rev. 25, 409–416. doi:10.1007/s10555-006-9005-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Bandettini, W. P., Kellman, P., Mancini, C., Booker, O. J., Vasu, S., Leung, S. W., et al. (2012). MultiContrast Delayed Enhancement (MCODE) Improves Detection of Subendocardial Myocardial Infarction by Late Gadolinium Enhancement Cardiovascular Magnetic Resonance: A Clinical Validation Study. J. Cardiovasc. Magn. Reson. 14, 1–10. doi:10.1186/1532-429X-14-83

PubMed Abstract | CrossRef Full Text | Google Scholar

Barrett, T., Wilhite, S. E., Ledoux, P., Evangelista, C., Kim, I. F., Tomashevsky, M., et al. (2013). NCBI GEO: Archive for Functional Genomics Data Sets-Update. Nucleic Acids Res. 41, D991–D995. doi:10.1093/nar/gks1193

PubMed Abstract | CrossRef Full Text | Google Scholar

Bhattacharya, S., Dunn, P., Thomas, C. G., Smith, B., Schaefer, H., Chen, J., et al. (2018). ImmPort, toward Repurposing of Open Access Immunological Assay Data for Translational and Clinical Research. Sci. Data 5, 1–9. doi:10.1038/sdata.2018.15

PubMed Abstract | CrossRef Full Text | Google Scholar

Brastianos, P. K., Shankar, G. M., Gill, C. M., Taylor-Weiner, A., Nayyar, N., Panka, D. J., et al. (2016). Dramatic Response of BRAF V600E Mutant Papillary Craniopharyngioma to Targeted Therapy. JNCI.J 108, djv310–5. doi:10.1093/jnci/djv310

PubMed Abstract | CrossRef Full Text | Google Scholar

Brastianos, P. K., Taylor-weiner, A., Manley, P. E., Jones, R. T., Dias-santagata, D., Thorner, A. R., et al. (2014). Exome Sequencing Identifies BRAF Mutations in Papillary Craniopharyngiomas. Nat. Genet. 46, 161–165. doi:10.1038/ng.2868.Exome

PubMed Abstract | CrossRef Full Text | Google Scholar

Bulk, E., Sargin, B., Krug, U., Hascher, A., Jun, Y., Knop, M., et al. (2009). S100A2 Induces Metastasis in Non-small Cell Lung Cancer. Clin. Cancer Res. 15, 22–29. doi:10.1158/1078-0432.CCR-08-0953

PubMed Abstract | CrossRef Full Text | Google Scholar

Cambien, B., Karimdjee, B. F., Richard-Fiardo, P., Bziouech, H., Barthel, R., Millet, M. A., et al. (2009). Organ-specific Inhibition of Metastatic colon Carcinoma by CXCR3 Antagonism. Br. J. Cancer 100, 1755–1764. doi:10.1038/sj.bjc.6605078

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, J., Tang, J., Chen, W., Gao, Y., He, Y., Zhang, Q., et al. (2017). Effects of Syndecan-1 on the Expression of Syntenin and the Migration of U251 Glioma Cells. Oncol. Lett. 14, 7217–7224. doi:10.3892/ol.2017.7170

PubMed Abstract | CrossRef Full Text | Google Scholar

Chin, C.-H., Chen, S.-H., Wu, H.-H., Ho, C.-W., Ko, M.-T., and Lin, C.-Y. (2014). cytoHubba: Identifying Hub Objects and Sub-networks from Complex Interactome. BMC Syst. Biol. 8, S11. doi:10.1186/1752-0509-8-S4-S11

PubMed Abstract | CrossRef Full Text | Google Scholar

Coffelt, S. B., Kersten, K., Doornebal, C. W., Weiden, J., Vrijland, K., Hau, C.-S., et al. (2015). IL-17-producing γδ T Cells and Neutrophils Conspire to Promote Breast Cancer Metastasis. Nature 522, 345–348. doi:10.1038/nature14282

PubMed Abstract | CrossRef Full Text | Google Scholar

Coy, S., Rashid, R., Lin, J.-R., Du, Z., Donson, A. M., Hankinson, T. C., et al. (2018). Multiplexed Immunofluorescence Reveals Potential PD-1/pd-L1 Pathway Vulnerabilities in Craniopharyngioma. Neuro. Oncol. 20, 1101–1112. doi:10.1093/neuonc/noy035

PubMed Abstract | CrossRef Full Text | Google Scholar

Daley, D., Zambirinisstantinos, C. P. P., Seifert, L., Akkad, N., Mohan, N., Werba, G., et al. (2016). Γδ T Cells Support Pancreatic Oncogenesis by Restraining αβ T Cell Activation. Cell 166, 1485–1499. e15. doi:10.1016/j.cell.2016.07.046

PubMed Abstract | CrossRef Full Text | Google Scholar

De Simone, V., Franzè, E., Ronchetti, G., Colantoni, A., Fantini, M. C., Di Fusco, D., et al. (2015). Th17-type Cytokines, IL-6 and TNF-α Synergistically Activate STAT3 and NF-kB to Promote Colorectal Cancer Cell Growth. Oncogene 34, 3493–3503. doi:10.1038/onc.2014.286

PubMed Abstract | CrossRef Full Text | Google Scholar

Dreyer, S. B., Pinese, M., Jamieson, N. B., Scarlett, C. J., Colvin, E. K., Pajic, M., et al. (2020). Precision Oncology in Surgery. Ann. Surg. 272, 366–376. doi:10.1097/SLA.0000000000003143

PubMed Abstract | CrossRef Full Text | Google Scholar

Elliott, R. E., Jane, J. A., and Wisoff, J. H. (2011). Surgical Management of Craniopharyngiomas in Children: Meta-Analysis and Comparison of Transcranial and Transsphenoidal Approaches. Neurosurgery 69, 630–643. doi:10.1227/NEU.0b013e31821a872d

PubMed Abstract | CrossRef Full Text | Google Scholar

Engebretsen, S., and Bohlin, J. (2019). Statistical Predictions with Glmnet. Clin. Epigenet 11, 10–12. doi:10.1186/s13148-019-0730-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Fleming, C., Morrissey, S., Cai, Y., and Yan, J. (2017). Γδ T Cells: Unexpected Regulators of Cancer Development and Progression. Trends Cancer 3, 561–570. doi:10.1016/j.trecan.2017.06.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Gong, J., Zhang, H., Xing, S., Li, C., Ma, Z., Jia, G., et al. (2014). High Expression Levels of CXCL12 and CXCR4 Predict Recurrence of Adamanti-Nomatous Craniopharyngiomas in Children. Cbm 14, 241–251. doi:10.3233/CBM-140397

PubMed Abstract | CrossRef Full Text | Google Scholar

Grob, S., Mirsky, D. M., Donson, A. M., Dahl, N., Foreman, N. K., Hoffman, L. M., et al. (2019). Targeting IL-6 Is a Potential Treatment for Primary Cystic Craniopharyngioma. Front. Oncol. 9, 1–6. doi:10.3389/fonc.2019.00791

PubMed Abstract | CrossRef Full Text | Google Scholar

Gump, J. M., Donson, A. M., Birks, D. K., Amani, V. M., Rao, K. K., Griesinger, A. M., et al. (2015). Identification of Targets for Rational Pharmacological Therapy in Childhood Craniopharyngioma. Acta Neuropathol. Commun. 3, 1–12. doi:10.1186/s40478-015-0211-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Hänzelmann, S., Castelo, R., and Guinney, J. (2013). GSVA: Gene Set Variation Analysis for Microarray and RNA-Seq Data. BMC Bioinformatics 14, 7. doi:10.1186/1471-2105-14-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Hegde, P. S., and Chen, D. S. (2020). Top 10 Challenges in Cancer Immunotherapy. Immunity 52, 17–35. doi:10.1016/j.immuni.2019.12.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, G., Zhang, J., Qing, G., Liu, D., Wang, X., Chen, Y., et al. (2021). S100A2 Silencing Relieves Epithelial-Mesenchymal Transition in Pulmonary Fibrosis by Inhibiting the Wnt/β-Catenin Signaling Pathway. DNA Cel Biol. 40, 18–25. doi:10.1089/dna.2020.6030

CrossRef Full Text | Google Scholar

Hussain, M., Adah, D., Tariq, M., Lu, Y., Zhang, J., and Liu, J. (2019). CXCL13/CXCR5 Signaling axis in Cancer. Life Sci. 227, 175–186. doi:10.1016/j.lfs.2019.04.053

PubMed Abstract | CrossRef Full Text | Google Scholar

Ito, K., and Murphy, D. (2013). Application of Ggplot2 to Pharmacometric Graphics. CPT: Pharmacometrics Syst. Pharmacol. 2, 79–16. doi:10.1038/psp.2013.56

PubMed Abstract | CrossRef Full Text | Google Scholar

Juratli, T. A., Jones, P. S., Wang, N., Subramanian, M., Aylwin, S. J. B., Odia, Y., et al. (2019). Targeted Treatment of Papillary Craniopharyngiomas Harboring BRAF V600E Mutations. Cancer 125, 2910–2914. doi:10.1002/cncr.32197

PubMed Abstract | CrossRef Full Text | Google Scholar

Kanehisa, M., Furumichi, M., Tanabe, M., Sato, Y., and Morishima, K. (2017). KEGG: New Perspectives on Genomes, Pathways, Diseases and Drugs. Nucleic Acids Res. 45, D353–D361. doi:10.1093/nar/gkw1092

PubMed Abstract | CrossRef Full Text | Google Scholar

Khafaga, Y., Jenkin, D., Kanaan, I., Hassounah, M., Al Shabanah, M., and Gray, A. (1998). Craniopharyngioma in Children. Int. J. Radiat. Oncology*Biology*Physics 42, 601–606. doi:10.1016/S0360-3016(98)00257-0

CrossRef Full Text | Google Scholar

Komohara, Y., Horlad, H., Ohnishi, K., Fujiwara, Y., Bai, B., Nakagawa, T., et al. (2012). Importance of Direct Macrophage - Tumor Cell Interaction on Progression of Human Glioma. Cancer Sci. 103, 2165–2172. doi:10.1111/cas.12015

PubMed Abstract | CrossRef Full Text | Google Scholar

Larkin, S. J., and Ansorge, O. (2013). Pathology and Pathogenesis of Craniopharyngiomas. Pituitary 16, 9–17. doi:10.1007/s11102-012-0418-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, J., Tang, Z., Wang, H., Wu, W., Zhou, F., Ke, H., et al. (2018). CXCL6 Promotes Non-small Cell Lung Cancer Cell Survival and Metastasis via Down-Regulation of miR-515-5p. Biomed. Pharmacother. 97, 1182–1188. doi:10.1016/j.biopha.2017.11.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, R., Liu, X., Zhou, X.-J., Chen, X., Li, J.-P., Yin, Y.-H., et al. (2020). Identification of a Prognostic Model Based on Immune-Related Genes of Lung Squamous Cell Carcinoma. Front. Oncol. 10, 1588. doi:10.3389/fonc.2020.01588

PubMed Abstract | CrossRef Full Text | Google Scholar

Liao, S., Liu, C., Zhu, G., Wang, K., Yang, Y., and Wang, C. (2020). Relationship between SDC1 and Cadherin Signalling Activation in Cancer. Pathol. - Res. Pract. 216, 152756. doi:10.1016/j.prp.2019.152756

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, D., Wang, Y., Zhou, Z., and Lin, Z. (2019). Immune Microenvironment of Primary and Recurrent Craniopharyngiomas: A Study of the Differences and Clinical Significance. World Neurosurg. 127, e212–e220. doi:10.1016/j.wneu.2019.02.212

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, J.-C., Sun, X.-W., Su, H., Chen, Q., Guo, T.-K., Li, Y., et al. (2017). Fibroblast-derived CXCL12/SDF-1α Promotes CXCL6 Secretion and Co-operatively Enhances Metastatic Potential through the PI3K/Akt/mTOR Pathway in colon Cancer. Wjg 23, 5167–5178. doi:10.3748/wjg.v23.i28.5167

PubMed Abstract | CrossRef Full Text | Google Scholar

Malgulwar, P. B., Nambirajan, A., Pathak, P., Faruq, M., Suri, V., Sarkar, C., et al. (2017). Study of β-catenin and BRAF Alterations in Adamantinomatous and Papillary Craniopharyngiomas: Mutation Analysis with Immunohistochemical Correlation in 54 Cases. J. Neurooncol. 133, 487–495. doi:10.1007/s11060-017-2465-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Mering, C. v., Huynen, M., Jaeggi, D., Schmidt, S., Bork, P., and Snel, B. (2003). STRING: A Database of Predicted Functional Associations between Proteins. Nucleic Acids Res. 31, 258–261. doi:10.1093/nar/gkg034

PubMed Abstract | CrossRef Full Text | Google Scholar

Mori, M., Takeshima, H., and Kuratsu, J.-i. (2004). Expression of Interleukin-6 in Human Craniopharyngiomas: a Possible Inducer of Tumor-Associated Inflammation. Int. J. Mol. Med. 14, 505–509. doi:10.3892/ijmm.14.4.505

PubMed Abstract | CrossRef Full Text | Google Scholar

Müller, H. L., Merchant, T. E., Warmuth-Metz, M., Martinez-Barbera, J.-P., and Puget, S. (2019). Craniopharyngioma. Nat. Rev. Dis. Primers 5, 75. doi:10.1038/s41572-019-0125-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Prabhu, V. C., and Brown, H. G. (2005). The Pathogenesis of Craniopharyngiomas. Childs Nerv Syst. 21, 622–627. doi:10.1007/s00381-005-1190-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Puchert, M., Obst, J., Koch, C., Zieger, K., and Engele, J. (2020). CXCL11 Promotes Tumor Progression by the Biased Use of the Chemokine Receptors CXCR3 and CXCR7. Cytokine 125, 154809. doi:10.1016/j.cyto.2019.154809

PubMed Abstract | CrossRef Full Text | Google Scholar

Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). Limma powers Differential Expression Analyses for RNA-Sequencing and Microarray Studies. Nucleic Acids Res. 43, e47. doi:10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

Roque, A., and Odia, Y. (2017). BRAF-V600E Mutant Papillary Craniopharyngioma Dramatically Responds to Combination BRAF and MEK Inhibitors. CNS Oncol. 6, 95–99. doi:10.2217/cns-2016-0034

PubMed Abstract | CrossRef Full Text | Google Scholar

Shankavaram, U. T., Varma, S., Kane, D., Sunshine, M., Chary, K. K., Reinhold, W. C., et al. (2009). CellMiner: A Relational Database and Query Tool for the NCI-60 Cancer Cell Lines. BMC Genomics 10, 277. doi:10.1186/1471-2164-10-277

PubMed Abstract | CrossRef Full Text | Google Scholar

Shannon, P., Markiel, A., Owen, O., Nitin, S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: A Software Environment for Integrated Models. Genome Res. 13, 2498–2504. doi:10.1101/gr.1239303.metabolite

PubMed Abstract | CrossRef Full Text | Google Scholar

Shi, S., Zhong, D., Xiao, Y., Wang, B., Wang, W., Zhang, F. a., et al. (2017a). Syndecan-1 Knockdown Inhibits Glioma Cell Proliferation and Invasion by Deregulating a C-src/FAK-Associated Signaling Pathway. Oncotarget 8, 40922–40934. doi:10.18632/oncotarget.16733

PubMed Abstract | CrossRef Full Text | Google Scholar

Shi, X. e., Zhou, Z., Wu, B., Zhang, Y., Qian, H., Sun, Y., et al. (2017b). Outcome of Radical Surgical Resection for Craniopharyngioma with Hypothalamic Preservation: A Single-Center Retrospective Study of 1054 Patients. World Neurosurg. 102, 167–180. doi:10.1016/j.wneu.2017.02.095

PubMed Abstract | CrossRef Full Text | Google Scholar

Tokunaga, R., Zhang, W., Naseem, M., Puccini, A., Berger, M. D., Soni, S., et al. (2018). CXCL9, CXCL10, CXCL11/CXCR3 axis for Immune Activation - a Target for Novel Cancer Therapy. Cancer Treat. Rev. 63, 40–47. doi:10.1016/j.ctrv.2017.11.007.CXCL9

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, H., Zhang, Z., Li, R., Ang, K. K., Zhang, H., Caraway, N. P., et al. (2005). Overexpression of S100A2 Protein as a Prognostic Marker for Patients with Stage I Non Small Cell Lung Cancer. Int. J. Cancer 116, 285–290. doi:10.1002/ijc.21035

PubMed Abstract | CrossRef Full Text | Google Scholar

Wightman, S. C., Uppal, A., Pitroda, S. P., Ganai, S., Burnette, B., Stack, M., et al. (2015). Oncogenic CXCL10 Signalling Drives Metastasis Development and Poor Clinical Outcome. Br. J. Cancer 113, 327–335. doi:10.1038/bjc.2015.193

PubMed Abstract | CrossRef Full Text | Google Scholar

Wolf, S., Haase-Kohn, C., and Pietzsch, J. (2011). S100A2 in Cancerogenesis: A Friend or a Foe? Amino Acids 41, 849–861. doi:10.1007/s00726-010-0623-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, G., Wang, L.-G., Han, Y., and He, Q.-Y. (2012). ClusterProfiler: An R Package for Comparing Biological Themes Among Gene Clusters. OMICS: A J. Integr. Biol. 16, 284–287. doi:10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

Zheng, S., Shen, T., Liu, Q., Liu, T., Tuerxun, A., Zhang, Q., et al. (2021). CXCL6 Fuels the Growth and Metastases of Esophageal Squamous Cell Carcinoma Cells Both In Vitro and In Vivo through Upregulation of PD‐L1 via Activation of STAT3 Pathway. J. Cel Physiol 236, 5373–5386. doi:10.1002/jcp.30236

CrossRef Full Text | Google Scholar