Research Paper Volume 15, Issue 24 pp 15535—15556
Potential impact of cuproptosis-related genes on tumor immunity in esophageal carcinoma
- 1 Department of Thoracic Surgery, The Second Hospital of Hebei Medical University, Shijiazhuang, China
- 2 Graduate school of Hebei Medical University, Shijiazhuang, China
- 3 Department of Thoracic Surgery, Hebei Chest Hospital, Shijiazhuang, China
- 4 Department of Thoracic Surgery, Affiliated Hospital of Chengde Medical University, Chengde, China
- 5 Clinical Laboratory Center, The Second Hospital of Hebei Medical University, Shijiazhuang, China
Received: July 10, 2023 Accepted: November 7, 2023 Published: December 30, 2023
https://doi.org/10.18632/aging.205391How to Cite
Copyright: © 2023 Guo et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Abstract
Cuproptosis involves a direct interaction with the tricarboxylic acid (TCA) lipid acylation components. This process intricately intersects with post-transcriptional lipid acylation (LA) and is linked to mitochondrial respiration and LA metabolism. Copper ions form direct bonds with acylated DLAT, promoting DLAT oligomerization, reducing Fe-S cluster proteins, and inducing a protein-triggered toxic stress response that culminates in cell demise. Simultaneously, the importance of immune contexture in cancer progression and treatment has significantly increased. We assessed the expression of cuproptosis-related genes (CRGs) across TCGA and validated our findings using the GEO data. Consensus clustering divided esophageal cancer (ESCA) patients into two clusters based on the expression of 7 CRGs. We evaluated the expression of immune checkpoint inhibitor (ICI) targets and calculated the elevated tumor mutational burden (TMB). Weighted gene co-expression network analysis (WGCNA) identified genes associated with the expression of CRGs and immunity. Cluster 1 exhibited increased immune infiltration, higher expression of ICI targets, higher TMB, and a higher incidence of deficiency in mismatch repair-microsatellite instability-high status. WGCNA analysis identified 14 genes associated with the expression of CRGs and immune scores. ROC analysis revealed specific hub genes with strong predictive capabilities. The expression levels of SLC6A3, MITD1, and PDHA1 varied across different pathological stages; CCS, LIPT2, PDHB, and PDHA1 showed variation in response to radiation therapy; MITD1 and PDHA1 exhibited differences related to the pathological M stages of ESCA. CRGs influence the immune contexture and can potentially transform cold tumors into hot tumors in ESCA patients.
Introduction
Esophageal carcinoma (ESCA) is globally recognized as the sixth most prevalent cancer, with a low survival rate [1]. Tumor recurrence and metastasis are identified as the primary causes of mortality among ESCA patients; however, the underlying mechanisms remain unclear [2]. ESCA exhibits significant variations in incidence, mortality, and histopathology across different geographic regions, especially in Southeast Asia/Africa and East Asia, contributing to a significant disease burden [3]. The two main histologic subtypes of ESCA are esophageal adenocarcinoma (EAC) and esophageal squamous cell carcinoma (ESCC), which together constitute over 95% of all ESCA cases [4]. Significant progress has been achieved in the treatment of ESCA in recent decades [5]. Combined approaches, including surgery, chemotherapy, and radiotherapy, have led to improved overall survival rates for patients with this disease [6]. However, a substantial number of ESCA patients are diagnosed at an advanced stage, frequently with delays in initial detection, resulting in a relatively poor prognosis and contributing to over 500,000 annual fatalities [7].
Copper, an essential micronutrient, plays a vital role in fundamental biological processes in all organisms. Being a redox-active metal, copper has the ability to donate and accept electrons, facilitating its transition between the reduced (Cu+) and oxidized (Cu2+) states [8]. Disruptions in copper homeostasis can lead to cellular toxicity, and changes in intracellular copper levels have been implicated in cancer development and progression [9]. Various copper ion carriers, including disulfide shimron, dithiocarbamate ester, chlorine compounds, and copper chelating agents like Sanchoishi and tetrathiomolybdenate, have been employed in cancer treatment using this mechanism [10–12]. Cuproptosis, a unique pathway, has become a significant player in the development of various cell death mechanisms, including apoptosis, pyroptosis, necroptosis, and ferroptosis [13]. The interaction between copper and the tricarboxylic acid (TCA) cycle influences the lipid acylation process, resulting in protein aggregation during acylation and the depletion of iron-sulfur cluster proteins. This leads to the accumulation of stress-inducing proteins and, ultimately, cell death [14]. Given copper’s dual function as both an essential enzyme cofactor and a potential inducer of cellular toxicity, it shows potential as an innovative therapeutic target. The accumulation of intracellular copper can be precisely targeted to eliminate cancer cells. Reports suggest that the combination of copper with platinum-based antitumor compounds can overcome drug resistance and function as a synergistic radiotherapeutic agent in cancer treatment [15].
Further research has unveiled a strong connection between cuproptosis and esophageal carcinoma. Irregular copper buildup is a common occurrence in different cancer types, which establishes a link between copper levels and cancer advancement [16, 17]. Copper plays a crucial role in promoting tumor growth and angiogenesis, acting as a cofactor for multiple pro-angiogenic molecules, including vascular endothelial growth factor (VEGF) [18]. Moreover, the copper ion plays an active role in the BRAF signaling pathway in cancer, promoting the proliferation and migration of tumor cells [19]. Notably, specific copper chelating agents and inhibitors, like ATOX1 and CCS, have shown the capacity to hinder the proliferation of various cancer cell types [18]. Tetrathiomolybdenate (TM), a copper chelating agent, has demonstrated good tolerability in both animal models and clinical trials, proving its effectiveness in inhibiting angiogenesis and tumor growth [20]. Therefore, the use of copper strips is considered a promising therapeutic approach for treating copper-rich cancers [21]. Nonetheless, the specific genes’ roles in cuproptosis and esophageal cancer have not been investigated. In this study, we pioneer the use of bioinformatics analysis to explore gene differences and assess the impact of genetic variations using gene enrichment analysis.
Results
Expression and consensus clustering of CRGs in ESCA
We performed a comparative analysis of the expression levels of 27 CRGs using a dataset that included 11 normal tissues and 152 tumor tissues from the TCGA database. For evaluating the proportions of stromal and immune cells in each sample, we computed four ESTIMATE indices. We aimed to investigate the role of CRGs in tumor immunity among ESCA patients. We analyzed the correlation between regulator expression and the ESTIMATE results (Figure 1A, 1B). Out of the identified regulators, six genes (SLC6A3, MITD1, CCS, LIPT2, ATOX1, and GLS) showed increased expression in tumor tissues, while PDHB had higher expression in normal tissues. The remaining 20 genes did not display significant differences in expression (Figure 1C). Next, we conducted a consensus clustering analysis on the expression matrix of CRGs using the 152 TCGA-ESCA samples. This analysis led to the classification of samples into two distinct clusters (Figure 1D and Supplementary Figures 1, 2). To explore the interactions among the 27 CRGs, we created a PPI network using online tools available on the STRING and GeneMANIA websites. The network analysis uncovered strong associations among the seven CRGs (Figure 1E, 1F).
Figure 1. Identification of CRGs and their association with the immune score and clustering of TCGA-ESCA. (A) Association between 27 CRGs and the ESTIMATE results. (B) Association between 7 CRGs and the ESTIMATE results. (C) Comparison of expression levels of the 7 CRGs in tumor and normal tissues. (D) TCGA-ESCA patients were classified into two clusters based on their expression levels of the 7 CRGs. (E) A PPI network of the differentially expressed genes associated with the 7 CRGs. (F) A PPI network of the differentially expressed genes associated with the 27 CRGs. *P < 0.05, **P < 0.01, ***P < 0.001.
Immunity and pathway enrichment analysis
We employed ESTIMATE, CIBERSORT, and ssGSEA analyses to understand the differences in immunological function between the two clusters. Results from the ESTIMATE algorithm showed that Cluster 1 had higher scores in stromal, immune, and estimate parameters and lower tumor purity compared to Cluster 2 (Figure 2A). Additionally, the CIBERSORT analysis revealed a higher proportion of CD4 T cells in Cluster 1 (Figure 2B). Conversely, ssGSEA analysis revealed increased expression of 25 immune cell subtypes, including immature B cells, macrophages, myeloid-derived suppressor cells (MDSCs), plasmacytoid dendritic cells, T follicular helper cells, and type 1 T helper cells, in Cluster 2 (Figure 2C). To explore the functional differences between the two clusters, we conducted GSEA with all the differentially expressed genes from Cluster 1 and Cluster 2. Importantly, we found several significant pathways related to immunity in the enrichment analysis using the MSigDB collection, one of which is the humoral immune response (Figure 2D).
Figure 2. Comparative analysis of immune characteristics between two clusters. (A) Expression level of stromal score, immune score, ESTIMATE score and tumor purity between the Cluster 1 and Cluster 2. (B) Comparison of immune cell proportions and (C) expression of immune cells between the two clusters. (D) Comparative analysis of functional enrichment. ns, not significant, *P < 0.05, **P < 0.01, ***P < 0.001.
Evaluation of sensitivity to immunotherapy
To assess the potential responsiveness of ESCA patients to immunotherapy, we examined the immunomodulatory drug targets identified in clinical trials for metastatic esophageal cancer. Then, we compared the expression levels of these immunomodulatory targets between the two clusters. Remarkably, we observed significantly higher expression of most targets, including PDCD1, CD274, PDCD1LG2, CTLA4, CD80, CD86, LAG3, HAVCR2, TIGHT, LGALS9, LAIR1, TNFRSF4, TNFRSF9, ICOS, CD40, and CD70, in Cluster 1 (Figure 3A–3D). These findings suggest that Cluster 1 may exhibit a more favorable response to immunotherapy than Cluster 2.
Figure 3. Comparative analysis of targets of immunomodulatory drugs in clinical trials for metastatic esophageal cancer between two clusters. The expression levels of immunomodulatory targets related to PD1 (A), CTLA4 (B), other immune checkpoint molecules (C), and agonists of T cell activation (D) varied between Cluster 1 and Cluster 2. *P < 0.05, **P < 0.01, ***P < 0.001.
Mutation analysis of CRGs
We acquired and analyzed mutation data from the TCGA database to study CRGs. Figure 4A shows that missense mutations are the most common variant. Among various types of variants, single nucleotide polymorphisms (SNPs) were the most common, with C > T being the predominant class of single nucleotide variant (SNV). TP53 notably showed the highest mutation frequency among the CRGs. Considering the substantial role of gene mutations in carcinogenesis, we delved into the distribution of somatic mutations in both Cluster 1 and Cluster 2. Figure 4B, 4C depict the corresponding results, showing that Cluster 1 had a somatic mutation rate of 98.78% (81 out of 82 samples), with missense mutations as the primary characteristic. Among these mutations, TP53 had the highest mutation frequency at 82%. In Cluster 2, the mutation rate was 100% (67 out of 67 samples), primarily consisting of missense mutations. TP53 exhibited the highest mutation frequency in this cluster, at 87%. TMB has become a potential marker for identifying cancer patients who may benefit from immunotherapy and predicting their therapeutic response to immune checkpoint inhibitors.
Figure 4. Comparative analysis of mutational landscapes between two clusters. (A) Overall mutational profile. Mutational landscape of Cluster 1 (B) and Cluster 2 (C).
WGCNA and identification of hub genes related with cuproptosis and immunity
We identified 987 differentially expressed genes, with 313 upregulated and 674 downregulated, between the two clusters. Their distribution was visualized using a volcano plot (Figure 5A). These genes were then subjected to WGCNA to investigate their co-expression patterns (Figure 5B, 5C). To identify a module associated with both cuproptosis and immunity, we conducted a correlation analysis between the modules and relevant traits (Figure 5D). From the available modules, we selected the black module due to its higher scores in stromal, immune, and estimate fractions, and lower tumor purity.
Figure 5. Identification of module genes associated with clustering and immunity in the WGCNA. (A) Volcano plot depicting differential analysis. (B) Analysis of network topology using soft powers. (C) Gene dendrogram with module colors. (D) Heatmap depicting the relationship between module eigengenes, clusters, and ESTIMATE results.
Functional enrichment of hub genes and their correlation with immune infiltration
To investigate the relationships among genes within the black module, we constructed a PPI network using the online tool provided on the STRING website. Figure 6A shows the strong associations among the 14 hub genes, which include CCL3L3, CCL5, CXCL11, CCL8, CXCL9, CXCL10, CXCL5, CCL7, CCL3, CXCL8, CD80, CSF3, CSF2, and FCGR2A. Additionally, we performed Spearman’s correlation analysis to investigate the relationships between these genes and immune infiltration, as evaluated using ESTIMATE and ssGSEA. The results indicated significant associations between the majority of genes and the immune response (Figure 6B, 6C).
Figure 6. Analysis of 14 hub genes. (A) PPI network of hub genes. (B) Correlation between hub genes and ESTIMATE results. (C) Correlation between hub genes and immune cell expression (ssGSEA).
GEO validation of immune characteristics between two clusters
Initially, we employed the same clustering approach to partition 42 esophageal cancer samples from the GSE199967 dataset into two clusters, mirroring the TCGA analysis (Figure 7A). We subsequently noted a comparable distribution of CRGs’ expression levels in these two clusters, which resembled the TCGA dataset’s findings. Subsequently, we evaluated the expression levels of immunomodulatory targets and the degree of immune infiltration using the CIBERSORT and ssGSEA methods (Figure 7B–7D). Cluster 1 notably exhibited significantly higher immune system activity in comparison to Cluster 2. Furthermore, we examined the expression of genes related to CRGs and immune checkpoint inhibitors in the GSE199967 cohort. We observed elevated expression levels of these genes in tumors compared to normal tissues (Figure 7E, 7F). We conducted a comprehensive evaluation of the potential implications of CRGs in immunotherapy and immune-related mechanisms. Our analysis showed that macrophages, including M0, M1, and M2 subsets, made up a significant portion of the GEO cohort (Figure 7G, 7H). In order to understand the mechanisms behind the various clinical risks associated with abnormal CRGs expression, we performed a GSEA analysis. The results revealed abnormal enrichment of pathways related to chromosomes and chromosome metabolism (Figure 7I).
Figure 7. Validation of immune contexture in GSE199967 between two clusters. (A) Division of GSE199967 patients into two clusters based on 7 CRGs. (B) Comparative analysis of targets of immunomodulatory drugs. (C) Proportion of immune cells and (D) expression of immune cells between the two clusters. (E) Comparative analysis of expression of 7 CRGs between tumor and normal tissues in GSE199967. (F) Comparative analysis of targets of immunomodulatory drugs between tumor and normal tissues in GSE199967. (G) Heatmap in GSE199967. (H) Percentages of immune cell types in GSE199967. (I) GSEA of the 7 CRGs. ns, not significant, *P < 0.05, **P < 0.01, ***P < 0.001.
External validation of 7 signature-associated genes
We performed qRT-PCR experiments to assess the expression levels of the 7 signature-associated genes in both esophageal cancer tissue and normal tissue. The results show a significant increase in the expression of SLC3A3, MITD1, LIPT2, ATOX1, PDHB, and GLS in tumor tissues compared to normal esophageal tissues (Figure 8A, 8B, 8D–8G). Of the six genes analyzed, PDHB showed an expression trend opposite to the decreasing trend observed in the TCGA database, while the other five genes exhibited a similar increasing trend as seen in the TCGA database. The expression level of CCS showed no significant difference between esophageal cancer and adjacent tissues (Figure 8C). The gene expression results from the clinical tissue samples closely matched the RNA sequencing data analyzed in the TCGA database.
Figure 8. Validation of expression of 7 CRGs in ESCA tissues using qRT-PCR. (A–G) The expression levels of the 7 CRGs in 10 pairs ESCA tissues and corresponding adjacent normal tissues were examined by qRT-PCR. ns, not significant, *P < 0.05, **P < 0.01, *** P < 0.001.
Prognosis prediction based on CRGs
We used a cohort of 151 patients with the TCGA-ESCA for external validation. Before conducting further analysis, we standardized the gene expression data using the “sva” package. Using relevant hub genes as predictive variables and taking the median risk score as a reference, we assessed the survival index of patients based on their outcome time and whether the outcome indicated death or survival (Figure 9A–9I). The ROC analysis showed the strong predictive ability of specific hub genes. Particularly, MITD1 displayed a positive predictive impact on the risk score, as supported by its AUC values of 0.6068 at 1 year, 0.5261 at 3 years, and 0.6727 at 5 years (Figure 9B). Similarly, PDHB and GLS exhibited strong predictive capabilities, with PDHB achieving AUC values of 0.6016 at 1 year, 0.5248 at 3 years, and 0.7788 at 5 years, and GLS achieving AUC values of 0.5627 at 1 year, 0.5943 at 3 years, and 0.7870 at 5 years (Figure 9F, 9G).
Figure 9. Evaluation of the independent prognostic value of gene expression using timeROC curves for 1-, 3-, and 5-year overall survival (OS) predictions of SLC6A3 (A), MITD1 (B), CCS (C), LIPT2 (D), ATOX1 (E), PDHB (F), GLS (G), DLAT (H), and DLD (I) through the nomogram in the TCGA cohort.
Differential expression of CRGs in different pathologic stages and histological grades of ESCA
As illustrated in Figure 10, the expression levels of three genes, namely SLC6A3, MITD1, and PDHA1, exhibited variations across different pathological stages (Figure 10A, 10B, 10G). Similarly, four genes, namely CCS, LIPT2, PDHB, and PDHA1, displayed variations in response to different radiation therapies (Figure 10D–10F, 10K). Of these genes, MITD1 and PDHA1 exhibited variations in different pathological M stages of ESCA. Specifically, the expression levels of SLC6A3, MITD1, and PDHA1 exhibited an increasing trend with tumor pathological stage, whereas CCS, LIPT2, PDHB, and PDHA1 displayed a decreasing trend in response to radiation therapy. However, MITD1 and PDHA1 displayed a contrasting trend in relation to the pathological M stage. Furthermore, the expression levels of PDHA1 varied across distinct pathological N stages, pathological T stages, and the count of metastatic lymph nodes in ESCA. Importantly, an increasing trend was observed in different pathological N stages and the count of metastatic lymph nodes, while a decreasing trend was observed in pathological T stages. These findings indicate a potential link between the expression of CRGs, disease grade, and radiotherapy resistance in ESCA. Patients with esophageal cancer were categorized into two groups, namely high expression and low expression, based on their PDHA1 expression. The Kaplan-Meier curve consistently showed that the high expression group had significantly shorter overall survival than the low expression group (Supplementary Figure 3).
Figure 10. Correlation between gene expression and clinicopathological staging characteristics. (A) Expression of SLC6A3 in different pathologic stages. (B, C) Expression of MITD1 in different pathologic stages and pathology M stage. (D) Expression of CCS, (E) LIPT2, and (F) PDHB in different radiation therapy. (G) Expression of PDHA1 in different pathologic stages, (H) pathology M stage, (I) pathology N stage, (J) pathology T stage, (K) radiation therapy, and (L) number of lymph nodes. *P < 0.05, **P < 0.01.
Discussion
This study is a pioneering investigation into the link between CRGs and ESCA. The main objective was to comprehensively analyze the expression patterns and immune infiltration of seven specific CRGs in patients with ESCA. Our analysis successfully identified two distinct molecular clusters based on the expression profiles of these CRGs. To clarify the functional differences between these clusters, we performed GSEA using all differentially expressed genes between Cluster 1 and Cluster 2. Our GSEA results unveiled several significant immunity-related pathways, including the humoral immune response. Importantly, we created a prognostic score that includes seven crucial CRGs (SLC6A3, MITD1, CCS, LIPT2, ATOX1, GLS, and PDHB). Of these genes, SLC6A3 has gained recognition as a risk factor due to its association with various familial mutations related to neuropsychiatric and neurological disorders [22]. Moreover, MITD1 deficiency has been discovered to hinder the growth and migration of clear cell renal cell carcinoma by inducing ferroptosis through the TAZ/SLC7A11 pathway [23]. Michael Grasso’s previous research showed that the copper chaperone for superoxide dismutase, CCS, binds temporarily to MEK1, promoting copper loading and thereby increasing MEK1/2 kinase activity [24]. The mechanism by which LIPT2 presence outside the mitochondrion induces apoptotic cell death is likely associated with mitochondrial dysfunction [25]. Inside human cells, ATOX1 plays a vital role in copper transport to the secretory pathway, where it transfers copper to copper-transporting ATPases (ATP7A and ATP7B) found in the trans-Golgi network and different endocytic vesicles. This mechanism promotes the maturation of copper-dependent enzymes in the secretory pathway and sustains copper levels in the cytosol and mitochondria [26]. Succinylation of GLS at K311 increases GLS activity, promoting glutaminolysis and generating NADPH and glutathione to counteract oxidative stress-induced reactive oxygen species (ROS) and apoptosis. As a result, this process facilitates tumor growth [27]. Manipulating PDHB expression through knockdown or overexpression resulted in decreased or increased expression of Myf5, MyoD, MyoG, and MyHC, resulting in reduced or increased rates of myogenic differentiation. Additionally, overexpressing PDHB in skeletal muscle alleviated D-galactose-induced sarcopenia in mice [28].
We assessed the immune characteristics of the two identified clusters using various methodologies, including ESTIMATE, CIBERSORT, and ssGSEA. ESTIMATE analysis showed that Cluster 1 had higher stromal, immune, and overall ESTIMATE scores, indicating a more dynamic tumor immune microenvironment. CIBERSORT analysis showed a significant increase in the proportion of CD4 memory activated T cells and resting NK cells in Cluster 1. Prior studies have emphasized that naive CD4 T lymphocytes differentiate into effector/memory cells during conventional adaptive immune responses upon recognizing foreign antigens [29]. In the immunosuppressive tumor microenvironment, NK cells may experience dysfunction due to exposure to inhibitory molecules produced by cancer cells, contributing to tumor escape [30]. Furthermore, the ssGSEA analysis indicated higher expression of 25 immune cell subtypes in Cluster 1, including CD8 T cells, CD4 T helper cells, dendritic cells (DCs), natural killer T (NKT) cells, regulatory T cells, and MDSCs. The infiltration of these immune cells has a significant impact on the clinical characteristics of ESCA. Higher intratumoral infiltration of CD8 T cells has been linked to longer survival times [31]. Th1 cells have been demonstrated to suppress ESCC cell proliferation, increase chemosensitivity and radiosensitivity, and correlate with improved prognosis [32]. DCs play a critical role in enhancing anti-tumor immunity by activating T cells [33]. Conventional type 1 dendritic cells are recruited into the tumor microenvironment upon being stimulated by NK cells [34]. Infiltration of LAMP-3-expressing DCs is positively correlated with intratumoral CD8 T cell levels and is linked to a favorable prognosis in ESCC [35]. NKT cells contribute to anti-tumor immunity by rejuvenating exhausted CD8 T cells in a tumor model resistant to anti-PD-1 therapy [36]. Macrophages are commonly categorized into M1 (proinflammatory; anti-tumor) and M2 (anti-inflammatory; tumor-promoting) subtypes [37]. Cluster 1 showed a higher proportion of M1 subtype macrophages, indicating a potential promotion of anti-tumor Th1-type responses, whereas Cluster 2 tended to establish a tolerogenic microenvironment. In summary, our analysis of the immune contexture showed increased immune cell infiltration in Cluster 1, suggesting higher immunological competence and the potential for immunotherapy benefits.
Previous studies extensively reported the FDA approval of ICIs that target programmed cell death 1 (PD-1), programmed cell death 1 ligand, and cytotoxic T lymphocyte antigen 4 [38]. Additionally, co-inhibitory receptor targets, like lymphocyte activation gene-3, T cell immunoglobulin-3, and T cell immunoglobulin and ITIM domain, have been identified as well [39]. This study aimed to compare two clusters of immunomodulatory drugs that have been studied in clinical trials for metastatic ESCC [40]. Importantly, Cluster 1 showed significantly higher expression levels of most of these targets. Subsequently, we analyzed the mutational profiles within the two clusters, revealing a significant difference between them. Specifically, Cluster 2 had a higher TMB compared to Cluster 1. TMB plays a critical role in influencing the generation of immunogenic peptides and, consequently, impacting the response to immunotherapy [41]. Therefore, these findings suggest that Cluster 1 may show a more favorable response to immunotherapy.
Subsequently, we conducted WGCNA to identify the blue module that exhibited distinct characteristics related to both CRGs and immune scores. By assessing module membership and gene significance values, we identified 14 genes within this module, which included CD80. CD80, among these genes, is associated with “professional antigen-presenting cells” that initiate T-cell activation via antigen presentation. Conversely, insufficient CD80 costimulation during antigen presentation may result in tolerance induction, and inhibiting CD80 costimulation has shown the potential to hinder the progression of autoimmune diseases in multiple animal models [42]. These occurrences are likely to be widespread in Cluster 1. Consequently, we hypothesize that Cluster 1 may demonstrate a more favorable response to immunotherapy than Cluster 2.
Cluster 1 exhibited “hot” tumor characteristics in its immunological features, while Cluster 2 showed characteristics indicative of a “cold” tumor [43]. By verifying CRGs expression in 10 pairs of cancer and normal tissues through qRT-PCR, we observed that Cluster 2 exhibited the overall immunological characteristics of ESCA, indicating a reduced immune response. It is plausible that CRGs may contribute to the conversion of a cold tumor to a hot tumor in ESCA. Furthermore, we developed a prognostic model based on CRGs and examined their expression in various pathological stages and histological grades of ESCA. Several hub genes showed robust predictive capabilities, as revealed by ROC analysis. Notably, PDHA1 expression displayed variations across distinct pathological N stages, T stages, and counts of metastatic lymph nodes in ESCA. Increasing trends were observed in N stages and lymph node count, whereas decreasing trends were noticed in T stages.
This study has contributed to our comprehension of the connection between CRGs and immunity; nonetheless, it’s important to recognize specific constraints. Firstly, the retrospective nature of the study underscores the importance of prioritizing future research on prospective studies, which can mitigate potential biases linked to retrospective designs. Furthermore, the lack of datasets with a sample size large enough and encompassing clinical prognostic information hampers the further validation of the results. Finally, the prognostic signature was built and confirmed using data from public databases. However, depending solely on TCGA and GEO databases restricts our capacity to explore CRGs’ expression at the protein level and reveal the direct mechanisms by which they influence anti-tumor immunity. Thus, it’s essential to collect more biological evidence to supplement the statistical results and enhance the robustness of the findings. Future research efforts should concentrate on unraveling the direct mechanisms that underlie the observed phenomena.
Materials and Methods
Data sources and preprocessing
Transcriptome profiling data, including HTSeq-Counts and HTSeq-FPKM, and clinical information, were acquired from the TCGA-ESCA project using R and the R package “TCGAbiolinks” [44]. The download was executed to collect extensive clinical data, which included age, sex, T stage, N stage, M stage, and prognostic information. Only cases with complete clinical information were included in the subsequent analyses. To further investigate, a log2(FPKM+1) transformation was applied to the Level 3 HTSeq-FPKM data of 173 primary solid tumor samples, while HTSeq-Counts were used for differential analysis.
To obtain nucleotide variation data, specifically MuTect2, for the 173 patients with ESCA, we used the R package “maftool” [45]. Mutational landscape analysis was exclusively conducted on these 173 patients due to the absence of mutation information for some ESCA patients. Waterfall plots, created with the R package “ComplexHeatmap” [46], were used to depict the genetic mutations observed in these patients. The tumor mutation burden (TMB) was calculated as the number of mutations per megabase using data on single nucleotide variations.
Expression profiling data from the GSE199967 array were obtained from the Gene Expression Omnibus (GEO) database. Subsequently, we extracted 27 cuproptosis-related genes (CRGs) from the reports by Tsvetkov [14] and the MsigDB. A total of 27 CRGs were included in the analysis conducted for this study (Supplementary Table 1). Seven specific CRGs were filtered based on the comparison of expression levels between normal tissues and ESCA samples.
Immune infiltration analysis
The ESTIMATE method, which analyzes gene expression patterns in tumor samples, quantified stromal and immune cell proportions in ESCA patients to evaluate the tumor microenvironment (TME). This assessment included the evaluation of stromal score (reflecting stromal content), immune score (indicating the extent of immune cell infiltration), ESTIMATE score (a composite measure of stroma and immune components), and tumor purity. These analyses were conducted using the R package “estimate” [47].
To estimate cell composition, we used the CIBERSORT computational method, which analyzes gene expression profiles. The deconvolution algorithm was applied to determine the relative abundance of 22 immune cell types in each ESCA patient [48]. The sum of the fractions of these 22 immune cell types in each sample equaled 1.
To assess immune cell type infiltration levels, we used the ssGSEA method from the R package “GSVA” [49]. This method allowed us to evaluate the expression profiles of 28 published gene sets specific to immune cells, giving insights into the infiltration levels of these 28 immune cell types [50].
Consensus clustering analysis of CRGs and GSEA analysis
We extracted the expression profiles of CRGs and conducted coherent clustering with the “ConsensusClusterPlus”| R package [51]. This analysis resulted in the samples being partitioned into two distinct clusters. Additionally, we identified consensus molecular subtypes (CMS) for each sample using the “CMScaller” R package [52]. CMS is a robust classification system for ESCA with distinct features: CMS1 (immune), CMS2 (canonical), CMS3 (metabolic), and CMS4 (mesenchymal) [53]. We employed a Sankey diagram to visualize the relationship between the two clusters and CMS.
Next, we performed gene set enrichment analysis (GSEA) on the differentially classified risk groups using the “ggplot2” and “clusterProfiler” packages in R software [54]. Pathway enrichment was deemed significant when it met specific criteria, including a normalized enrichment score (|NES| > 1), a P-value < 0.05, and a false discovery rate (FDR) q-value < 0.05.
Differential expressed genes
Expression profiling data (HTSeq-Counts) were analyzed with the “DESeq2” R package [55] to identify differentially expressed genes (DEGs) between the two clusters. The criteria for selecting DEGs included a threshold of |log2FoldChange| > 1 and an adjusted P-value < 0.05.
Weighted gene co-expression network analysis and functional enrichment analysis
Weighted gene co-expression network analysis (WGCNA) was performed on the differentially expressed genes with the “WGCNA” R package [56]. To ensure the construction of a co-expression network with a scale-free distribution, we chose a soft power of 5. Eleven modules were identified, and we assessed their associations with the cluster, stromal score, immune score, ESTIMATE score, and tumor purity. Subsequently, we identified 14 genes based on calculations of module membership (MM) and gene significance (GS).
We performed Gene Ontology (GO) analysis using the “clusterProfiler” R package [54] to gain insights into the functions of the 14 selected DEGs. Additionally, we constructed a protein-protein interaction (PPI) network using the STRING database [57]. Furthermore, we used the “corrplot” R package to conduct Spearman’s correlation analysis, investigating the relationships between genes, gene-ESTIMATE, and gene-ssGSEA.
Tissue samples and qRT-PCR
We obtained 10 tumor tissue samples and corresponding adjacent normal esophageal tissue samples from patients who underwent tumor resection for ESCA. The tissue sample collection took place at the Thoracic Surgery Department of the Second Hospital of Hebei Medical University, with approval from the hospital’s Medical Ethics Committee. To maintain sample integrity, fresh tumor and non-tumor tissues were promptly frozen in liquid nitrogen. RNA extraction utilized the TRIzol Reagent (Invitrogen, USA). cDNA synthesis was conducted using the PrimeScriptTM RT reagent Kit with gDNA Eraser (Takara, Beijing, China) via reverse transcription. Quantitative real-time polymerase chain reaction (qRT-PCR) analysis utilized the SYBR Premix Ex Taq (Takara). Expression data were normalized to the internal control GAPDH using the 2−ΔΔCT method. Gene-specific primers were synthesized by Sangon Biotech (Shanghai, China), and their sequences are available in Supplementary Table 2.
Prognosis prediction of CRGs and difference analysis of scores with clinical stages
Using the median risk score as the basis, we employed key hub genes as predictive variables to estimate patients’ survival index, taking into account their survival time and outcome (death or survival). Furthermore, we used the signature to calculate the 1-year, 3-year, and 5-year survival rates through the nearest neighbor estimation method [58]. To evaluate its predictive performance, we generated receiver operating characteristic (ROC) curves with the “survivalROC” R package. The clinicopathological data of ESCA samples were retrieved from the TCGA database. Statistical analysis was conducted using the R programming language. The significance was assessed by employing either the Wilcoxon rank sum test or the Kruskal-Wallis rank sum test for comparisons, depending on the number of clinical stages.
Statistical analysis
All statistical analyses were performed using R software version 4.2.1. Figures were created and compiled using Adobe Illustrator. The Wilcoxon rank-sum test was used for box plot analyses. Spearman’s coefficient was used for correlation analysis. The chi-square test was used to compare clinical characteristics between the two clusters, and Fisher’s exact test was used when necessary. Multivariate logistic regression analysis was performed to assess the impact of clinical characteristics on the clusters. All hypothesis tests were two-sided, and statistical significance was defined as a P-value less than 0.05.
Abbreviations
LA: lipid acylation; ESCA: esophageal carcinoma; EAC: esophageal adenocarcinoma; TCA: tricarboxylic acid; VEGF: vascular endothelial growth factor; BRAF: v-raf murine sarcoma viral oncogene homolog B1; TM: tetrathiomolybdenate; CRGs: cuproptosis-related genes; TCGA: The Cancer Genome Atlas; PPI: protein-protein interactions; MDSCs: myeloid-derived suppressor cells; GSEA: gene set enrichment analysis; SNPs: single nucleotide polymorphisms; SNV: single nucleotide variant; TMB: tumor mutational burden; WGCNA: weighted gene co-expression network analysis; GEO: Gene Expression Omnibus; AUC: area under curve; ROS: reactive oxygen species; DCs: dendritic cells; ICIs: immune checkpoint inhibitors; PD-1: targeting programmed cell death 1; TME: tumor microenvironment; GSVA: gene set variation analysis; DEGs: differentially expressed genes; GO: gene ontology.
Author Contributions
SL had complete access to the study data and assumed responsibility for data integrity and accuracy of the analysis. The concept and design were developed by SL. PG and ZN performed the experiment and drafted the manuscript. The manuscript was revised by SL. Statistical analysis was performed by PG and FZ. SL and ZL provided supervision. The final manuscript was read and approved by all authors.
Acknowledgments
The authors thank the Department of Thoracic Surgery and the Department of urology Surgery, The Second Hospital of Hebei Medical University and the contributions from the TCGA and GEO network.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Ethical Statement and Consent
The studies involving human participants were reviewed and approved by The Ethics Committee of the Second Hospital of Hebei Medical University (approval number 2022-R080). The patients/participants provided their written informed consent to participate in this study.
Funding
This study was supported by Hebei Provincial 2024 Government-Funded Clinical Medicine Excellent Talents Training Projects (No: ZF2024177 to Yishuai Li and 303-2022-27-37 to Shujun Li) and Hebei Province Science and Technology Support Program (No.223777106D to Zhirong Li).
References
- 1. Poston GJ. Global cancer surgery: The Lancet Oncology review. Eur J Surg Oncol. 2015; 41:1559–61. https://doi.org/10.1016/j.ejso.2015.09.004 [PubMed]
- 2. Hou X, Wei JC, Fu JH, Wang X, Luo RZ, He JH, Zhang LJ, Lin P, Yang HX. Vascular Endothelial Growth Factor is a Useful Predictor of Postoperative Distant Metastasis and Survival Prognosis in Esophageal Squamous Cell Carcinoma. Ann Surg Oncol. 2015; 22:3666–73. https://doi.org/10.1245/s10434-015-4390-x [PubMed]
- 3. Malhotra GK, Yanala U, Ravipati A, Follet M, Vijayakumar M, Are C. Global trends in esophageal cancer. J Surg Oncol. 2017; 115:564–79. https://doi.org/10.1002/jso.24592 [PubMed]
- 4. Abbas G, Krasna M. Overview of esophageal cancer. Ann Cardiothorac Surg. 2017; 6:131–6. https://doi.org/10.21037/acs.2017.03.03 [PubMed]
- 5. Caruz A, Samsom M, Alonso JM, Alcami J, Baleux F, Virelizier JL, Parmentier M, Arenzana-Seisdedos F. Genomic organization and promoter characterization of human CXCR4 gene. FEBS Lett. 1998; 426:271–8. https://doi.org/10.1016/s0014-5793(98)00359-7 [PubMed]
- 6. Wang B, Liu G, Ding L, Zhao J, Lu Y. FOXA2 promotes the proliferation, migration and invasion, and epithelial mesenchymal transition in colon cancer. Exp Ther Med. 2018; 16:133–40. https://doi.org/10.3892/etm.2018.6157 [PubMed]
- 7. Wegner SA, Ehrenberg PK, Chang G, Dayhoff DE, Sleeker AL, Michael NL. Genomic organization and functional characterization of the chemokine receptor CXCR4, a major entry co-receptor for human immunodeficiency virus type 1. J Biol Chem. 1998; 273:4754–60. https://doi.org/10.1074/jbc.273.8.4754 [PubMed]
- 8. Arredondo M, Núñez MT. Iron and copper metabolism. Mol Aspects Med. 2005; 26:313–27. https://doi.org/10.1016/j.mam.2005.07.010 [PubMed]
- 9. Babak MV, Ahn D. Modulation of Intracellular Copper Levels as the Mechanism of Action of Anticancer Copper Complexes: Clinical Relevance. Biomedicines. 2021; 9:852. https://doi.org/10.3390/biomedicines9080852 [PubMed]
- 10. Brady DC, Crowe MS, Greenberg DN, Counter CM. Copper Chelation Inhibits BRAFV600E-Driven Melanomagenesis and Counters Resistance to BRAFV600E and MEK1/2 Inhibitors. Cancer Res. 2017; 77:6240–52. https://doi.org/10.1158/0008-5472.CAN-16-1190 [PubMed]
- 11. Davis CI, Gu X, Kiefer RM, Ralle M, Gade TP, Brady DC. Altered copper homeostasis underlies sensitivity of hepatocellular carcinoma to copper chelation. Metallomics. 2020; 12:1995–2008. https://doi.org/10.1039/d0mt00156b [PubMed]
- 12. Chen D, Cui QC, Yang H, Dou QP. Disulfiram, a clinically used anti-alcoholism drug and copper-binding agent, induces apoptotic cell death in breast cancer cultures and xenografts via inhibition of the proteasome activity. Cancer Res. 2006; 66:10425–33. https://doi.org/10.1158/0008-5472.CAN-06-2126 [PubMed]
- 13. Huang H, Chen Z, Sun Q. Mammalian Cell Competitions, Cell-in-Cell Phenomena and Their Biomedical Implications. Curr Mol Med. 2015; 15:852–60. https://doi.org/10.2174/1566524015666151026101101 [PubMed]
- 14. Tsvetkov P, Coy S, Petrova B, Dreishpoon M, Verma A, Abdusamad M, Rossen J, Joesch-Cohen L, Humeidi R, Spangler RD, Eaton JK, Frenkel E, Kocak M, et al. Copper induces cell death by targeting lipoylated TCA cycle proteins. Science. 2022; 375:1254–61. https://doi.org/10.1126/science.abf0529 [PubMed]
- 15. Ge EJ, Bush AI, Casini A, Cobine PA, Cross JR, DeNicola GM, Dou QP, Franz KJ, Gohil VM, Gupta S, Kaler SG, Lutsenko S, Mittal V, et al. Connecting copper and cancer: from transition metal signalling to metalloplasia. Nat Rev Cancer. 2022; 22:102–13. https://doi.org/10.1038/s41568-021-00417-2 [PubMed]
- 16. Margalioth EJ, Schenker JG, Chevion M. Copper and zinc levels in normal and malignant tissues. Cancer. 1983; 52:868–72. https://doi.org/10.1002/1097-0142(19830901)52:5<868::aid-cncr2820520521>3.0.co;2-k [PubMed]
- 17. Ishida S, McCormick F, Smith-McCune K, Hanahan D. Enhancing tumor-specific uptake of the anticancer drug cisplatin with a copper chelator. Cancer Cell. 2010; 17:574–83. https://doi.org/10.1016/j.ccr.2010.04.011 [PubMed]
- 18. Barresi V, Trovato-Salinaro A, Spampinato G, Musso N, Castorina S, Rizzarelli E, Condorelli DF. Transcriptome analysis of copper homeostasis genes reveals coordinated upregulation of SLC31A1, SCO1, and COX11 in colorectal cancer. FEBS Open Bio. 2016; 6:794–806. https://doi.org/10.1002/2211-5463.12060 [PubMed]
- 19. Finney L, Vogt S, Fukai T, Glesne D. Copper and angiogenesis: unravelling a relationship key to cancer progression. Clin Exp Pharmacol Physiol. 2009; 36:88–94. https://doi.org/10.1111/j.1440-1681.2008.04969.x [PubMed]
- 20. Wang J, Luo C, Shan C, You Q, Lu J, Elf S, Zhou Y, Wen Y, Vinkenborg JL, Fan J, Kang H, Lin R, Han D, et al. Inhibition of human copper trafficking by a small molecule significantly attenuates cancer cell proliferation. Nat Chem. 2015; 7:968–79. https://doi.org/10.1038/nchem.2381 [PubMed]
- 21. Chan N, Willis A, Kornhauser N, Ward MM, Lee SB, Nackos E, Seo BR, Chuang E, Cigler T, Moore A, Donovan D, Vallee Cobham M, Fitzpatrick V, et al. Influencing the Tumor Microenvironment: A Phase II Study of Copper Depletion Using Tetrathiomolybdate in Patients with Breast Cancer at High Risk for Recurrence and in Preclinical Models of Lung Metastases. Clin Cancer Res. 2017; 23:666–76. https://doi.org/10.1158/1078-0432.CCR-16-1326 [PubMed]
- 22. Reith MEA, Kortagere S, Wiers CE, Sun H, Kurian MA, Galli A, Volkow ND, Lin Z. The dopamine transporter gene SLC6A3: multidisease risks. Mol Psychiatry. 2022; 27:1031–46. https://doi.org/10.1038/s41380-021-01341-5 [PubMed]
- 23. Zhang Y, Li Y, Qiu Q, Chen Z, Du Y, Liu X. MITD1 Deficiency Suppresses Clear Cell Renal Cell Carcinoma Growth and Migration by Inducing Ferroptosis through the TAZ/SLC7A11 Pathway. Oxid Med Cell Longev. 2022; 2022:7560569. https://doi.org/10.1155/2022/7560569 [PubMed]
- 24. Grasso M, Bond GJ, Kim YJ, Boyd S, Matson Dzebo M, Valenzuela S, Tsang T, Schibrowsky NA, Alwan KB, Blackburn NJ, Burslem GM, Wittung-Stafshede P, Winkler DD, et al. The copper chaperone CCS facilitates copper binding to MEK1/2 to promote kinase activation. J Biol Chem. 2021; 297:101314. https://doi.org/10.1016/j.jbc.2021.101314 [PubMed]
- 25. Bernardinelli E, Costa R, Scantamburlo G, To J, Morabito R, Nofziger C, Doerrier C, Krumschnabel G, Paulmichl M, Dossena S. Mis-targeting of the mitochondrial protein LIPT2 leads to apoptotic cell death. PLoS One. 2017; 12:e0179591. https://doi.org/10.1371/journal.pone.0179591 [PubMed]
- 26. Hatori Y, Lutsenko S. The Role of Copper Chaperone Atox1 in Coupling Redox Homeostasis to Intracellular Copper Distribution. Antioxidants (Basel). 2016; 5:25. https://doi.org/10.3390/antiox5030025 [PubMed]
- 27. Tong Y, Guo D, Lin SH, Liang J, Yang D, Ma C, Shao F, Li M, Yu Q, Jiang Y, Li L, Fang J, Yu R, Lu Z. SUCLA2-coupled regulation of GLS succinylation and activity counteracts oxidative stress in tumor cells. Mol Cell. 2021; 81:2303–16.e8. https://doi.org/10.1016/j.molcel.2021.04.002 [PubMed]
- 28. Jiang X, Ji S, Yuan F, Li T, Cui S, Wang W, Ye X, Wang R, Chen Y, Zhu S. Pyruvate dehydrogenase B regulates myogenic differentiation via the FoxP1-Arih2 axis. J Cachexia Sarcopenia Muscle. 2023; 14:606–21. https://doi.org/10.1002/jcsm.13166 [PubMed]
- 29. Kawabe T, Sher A. Memory-phenotype CD4+ T cells: a naturally arising T lymphocyte population possessing innate immune function. Int Immunol. 2022; 34:189–96. https://doi.org/10.1093/intimm/dxab108 [PubMed]
- 30. Cózar B, Greppi M, Carpentier S, Narni-Mancinelli E, Chiossone L, Vivier E. Tumor-Infiltrating Natural Killer Cells. Cancer Discov. 2021; 11:34–44. https://doi.org/10.1158/2159-8290.CD-20-0655 [PubMed]
- 31. Liu B, Liu Z, Gao C. Relationship Between CD8+ T Cells and Prognosis of Esophageal Cancer Patients: A Systematic Review and Meta-analysis. Mol Biotechnol. 2023. [Epub ahead of print]. https://doi.org/10.1007/s12033-023-00733-y [PubMed]
- 32. Yuan J, Weng Z, Tan Z, Luo K, Zhong J, Xie X, Qu C, Lin X, Yang H, Wen J, Fu J. Th1-involved immune infiltrates improve neoadjuvant chemoradiotherapy response of esophageal squamous cell carcinoma. Cancer Lett. 2023; 553:215959. https://doi.org/10.1016/j.canlet.2022.215959 [PubMed]
- 33. Wculek SK, Cueto FJ, Mujal AM, Melero I, Krummel MF, Sancho D. Dendritic cells in cancer immunology and immunotherapy. Nat Rev Immunol. 2020; 20:7–24. https://doi.org/10.1038/s41577-019-0210-z [PubMed]
- 34. Böttcher JP, Bonavita E, Chakravarty P, Blees H, Cabeza-Cabrerizo M, Sammicheli S, Rogers NC, Sahai E, Zelenay S, Reis e Sousa C. NK Cells Stimulate Recruitment of cDC1 into the Tumor Microenvironment Promoting Cancer Immune Control. Cell. 2018; 172:1022–37.e14. https://doi.org/10.1016/j.cell.2018.01.004 [PubMed]
- 35. Nishimura J, Tanaka H, Yamakoshi Y, Hiramatsu S, Tamura T, Toyokawa T, Muguruma K, Maeda K, Hirakawa K, Ohira M. Impact of tumor-infiltrating LAMP-3 dendritic cells on the prognosis of esophageal squamous cell carcinoma. Esophagus. 2019; 16:333–44. https://doi.org/10.1007/s10388-019-00669-w [PubMed]
- 36. Bae EA, Seo H, Kim BS, Choi J, Jeon I, Shin KS, Koh CH, Song B, Kim IK, Min BS, Han YD, Shin SJ, Kang CY. Activation of NKT Cells in an Anti-PD-1-Resistant Tumor Model Enhances Antitumor Immunity by Reinvigorating Exhausted CD8 T Cells. Cancer Res. 2018; 78:5315–26. https://doi.org/10.1158/0008-5472.CAN-18-0734 [PubMed]
- 37. Aras S, Zaidi MR. TAMeless traitors: macrophages in cancer progression and metastasis. Br J Cancer. 2017; 117:1583–91. https://doi.org/10.1038/bjc.2017.356 [PubMed]
-
38.
Ganesh K, Stadler ZK, Cercek A, Mendelsohn RB, Shia J, Segal NH, Diaz LA
Jr . Immunotherapy in colorectal cancer: rationale, challenges and potential. Nat Rev Gastroenterol Hepatol. 2019; 16:361–75. https://doi.org/10.1038/s41575-019-0126-x [PubMed] - 39. Anderson AC, Joller N, Kuchroo VK. Lag-3, Tim-3, and TIGIT: Co-inhibitory Receptors with Specialized Functions in Immune Regulation. Immunity. 2016; 44:989–1004. https://doi.org/10.1016/j.immuni.2016.05.001 [PubMed]
- 40. Thuss-Patience P, Stein A. Immunotherapy in Squamous Cell Cancer of the Esophagus. Curr Oncol. 2022; 29:2461–71. https://doi.org/10.3390/curroncol29040200 [PubMed]
- 41. Havel JJ, Chowell D, Chan TA. The evolving landscape of biomarkers for checkpoint inhibitor immunotherapy. Nat Rev Cancer. 2019; 19:133–50. https://doi.org/10.1038/s41568-019-0116-x [PubMed]
- 42. Lu P, Wang YL, Linsley PS. Regulation of self-tolerance by CD80/CD86 interactions. Curr Opin Immunol. 1997; 9:858–62. https://doi.org/10.1016/s0952-7915(97)80190-2 [PubMed]
- 43. Galon J, Bruni D. Approaches to treat immune hot, altered and cold tumours with combination immunotherapies. Nat Rev Drug Discov. 2019; 18:197–218. https://doi.org/10.1038/s41573-018-0007-y [PubMed]
- 44. Colaprico A, Silva TC, Olsen C, Garofano L, Cava C, Garolini D, Sabedot TS, Malta TM, Pagnotta SM, Castiglioni I, Ceccarelli M, Bontempi G, Noushmehr H. TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res. 2016; 44:e71. https://doi.org/10.1093/nar/gkv1507 [PubMed]
- 45. Mayakonda A, Lin DC, Assenov Y, Plass C, Koeffler HP. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res. 2018; 28:1747–56. https://doi.org/10.1101/gr.239244.118 [PubMed]
- 46. Gu Z, Eils R, Schlesner M. Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics. 2016; 32:2847–9. https://doi.org/10.1093/bioinformatics/btw313 [PubMed]
- 47. Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, Treviño V, Shen H, Laird PW, Levine DA, Carter SL, Getz G, Stemke-Hale K, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013; 4:2612. https://doi.org/10.1038/ncomms3612 [PubMed]
- 48. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, Hoang CD, Diehn M, Alizadeh AA. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015; 12:453–7. https://doi.org/10.1038/nmeth.3337 [PubMed]
- 49. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013; 14:7. https://doi.org/10.1186/1471-2105-14-7 [PubMed]
- 50. Bindea G, Mlecnik B, Tosolini M, Kirilovsky A, Waldner M, Obenauf AC, Angell H, Fredriksen T, Lafontaine L, Berger A, Bruneval P, Fridman WH, Becker C, et al. Spatiotemporal dynamics of intratumoral immune cells reveal the immune landscape in human cancer. Immunity. 2013; 39:782–95. https://doi.org/10.1016/j.immuni.2013.10.003 [PubMed]
- 51. Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics. 2010; 26:1572–3. https://doi.org/10.1093/bioinformatics/btq170 [PubMed]
- 52. Eide PW, Bruun J, Lothe RA, Sveen A. CMScaller: an R package for consensus molecular subtyping of colorectal cancer pre-clinical models. Sci Rep. 2017; 7:16618. https://doi.org/10.1038/s41598-017-16747-x [PubMed]
- 53. Guinney J, Dienstmann R, Wang X, de Reyniès A, Schlicker A, Soneson C, Marisa L, Roepman P, Nyamundanda G, Angelino P, Bot BM, Morris JS, Simon IM, et al. The consensus molecular subtypes of colorectal cancer. Nat Med. 2015; 21:1350–6. https://doi.org/10.1038/nm.3967 [PubMed]
- 54. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012; 16:284–7. https://doi.org/10.1089/omi.2011.0118 [PubMed]
- 55. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014; 15:550. https://doi.org/10.1186/s13059-014-0550-8 [PubMed]
- 56. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008; 9:559. https://doi.org/10.1186/1471-2105-9-559 [PubMed]
- 57. Szklarczyk D, Gable AL, Lyon D, Junge A, Wyder S, Huerta-Cepas J, Simonovic M, Doncheva NT, Morris JH, Bork P, Jensen LJ, Mering CV. STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res. 2019; 47:D607–13. https://doi.org/10.1093/nar/gky1131 [PubMed]
- 58. Heagerty PJ, Zheng Y. Survival model predictive accuracy and ROC curves. Biometrics. 2005; 61:92–105. https://doi.org/10.1111/j.0006-341X.2005.030814.x [PubMed]