Comprehensive Analysis of TRP Channel-Related Genes for Estimating the Immune Microenvironment, Prognosis, and Therapeutic Effect in Patients With Esophageal Squamous Cell Carcinoma

RNA sequencing data and clinical information of ESCC patients were obtained from the TCGA database (80 tumor samples and 11 normal samples) and the GEO database (GSE35624↗, including 119 tumor samples). Meanwhile, the copy number variation (CNV), single nucleotide variation (SNV) and methylation data of ESCC were downloaded from TCGA database. To collect TRP channel-related genes, we downloaded the Reactome_TRP_channels gene set from the MSigDB database, and another gene set from the KEGG database: inflammatory mediator regulation of TRP channels. Finally, we obtained 120 TRP channel-related genes for subsequent analyses.

Differential expression of TRP channel-related genes in tumor and normal tissues was detected using the “limma” R package. The “maftools” package was used to generate a waterfall plot of the mutation frequency of differentially expressed genes (DEGs) in ESCC patients. A protein-protein interaction (PPI) network was mapped with online tools on the GeneMANIA website. In addition, the gene set variation analysis (GSVA) of R software (Hänzelmann et al., 2013) was employed to predict the activation or inhibition status of TRP channel-related DEGs and cancer-associated pathways.

Two cohorts were used to validate the prediction performance while one of the cohorts (GEO-ESCC, GSE35624↗) was used to construct the prognostic signature. Significant prognostic genes were identified by univariate Cox regression analysis (p < 0.1). Subsequently, we performed the least absolute shrinkage and selection operator (LASSO) Cox regression analysis and multivariate Cox regression analysis using the “glmnet” package to screen for genes involved in the TRP scoring formula. The TRP risk score was calculated as: (gene 1 × expression coefficient) + (gene 2 × expression coefficient) + … + (gene n × expression coefficient). All patients in TCGA and GEO were divided into two groups (low-risk group and high-risk groups) according to the median value of TRP risk score. The overall survival (OS) between the high-risk and low-risk groups was compared with Kaplan-Meier analysis. In addition, principal component analysis (PCA) analysis was used to determine the ability of the TRP risk score to discriminate between patients. The receiver operating characteristic (ROC) analysis was performed using the “timeROC” R package. Cox regression, both univariate and multivariate, were used to identify independent prognostic factors in patients with ESCC. The nomogram was created and their performance was evaluated using the “rms” R package.

The CIBERSORT algorithm calculates the proportion of different immune cell types based on the expression levels of immune cell-related genes. The output of the 22 infiltrated immune cells was integrated to generate an immune cell matrix for subsequent analysis. We performed Gene set enrichment analyses (GSEA) to identify 672 GO terms and 30 KEGG pathways associated with different risk groups. In addition, we used the “pRRophetic” package to predict the sensitivity of different risk groups to chemotherapy.

A total of 10 tumor tissue samples and nearby normal esophageal tissue samples were obtained from ESCC patients who underwent tumor resection. All tissue samples were collected from the Thoracic Surgery Department of the Second Hospital of Hebei Medical University with the approval by the Medical Ethics Committee of the hospital. In liquid nitrogen, fresh tumor and non-tumor tissues were frozen in a snap, preserving their integrity. The extraction of RNA was carried out using the TRIzol Reagent (Invitrogen, United States). PrimeScriptTM RT reagent Kit with gDNA Eraser (Takara) was adopted to synthesize complementary DNA (cDNA) by reverse transcription reaction. The SYBR Premix Ex Taq (Takara) was adopted for the qRT-PCR analysis. The normalization of all expression data to GAPDH as an internal control was made with the 2^−ΔΔCT approach. Sangon Biotech (Sangon Biotech, Shanghai, China) was applied to chemically synthesize all primers. Supplementary File S1 lists the primer sequences and original data. The score was further standardized and simplified to generate a riskscore. The score was subsequently mapped by subtracting the minimum and dividing by the maximum. Mapping was conducted to facilitate the interpretation of results from different platforms. The riskscore was calculated as follows: Riskscore (qRT-PCR) =(Score-Min)/Max.

The R software (v.4.0.1) was employed to perform all statistical analyses. The bioinformatics approach section covers detailed statistical approaches about Transcriptome information. p < 0.05 was of statistical significance.

The flow chart of this study is shown in Supplementary Figure S1. We obtained a gene set (Reactom_TRP_channels) from the MSigDB database, and another gene set (inflammatory mediator regulation of TRP channels) from the KEGG database, as shown in Figure 1A. Finally, 120 TRP channel-related genes were collected for subsequent analyses. RNA-seq data from the TCGA-ESCC cohort were analyzed using the limma package (p < 0.05; |logFC| > 1.5), and a total of 18 TRP channel-associated DEGs were identified in 80 ESCC and 11 normal samples (Figure 1B). In the heat map, the genes marked in red were members of the TRP superfamily, including TRPM3, TRPM35, TRPV4, and TRPV3. In addition, the volcano map showed that 10 genes were upregulated and 8 genes were downregulated in tumor samples compared to normal samples (Figure 1C). Specifically, the 10 upregulated genes were: TRPV4, TRPV3, HTR2C, PLA2G4E, ASIC5, CALML3, CALML5, IL1RAP, MAPK12, and NGF. 8 downregulated genes were TRPM5, TRPM3, PLCB4, ITPR1, PTGER2, ADCY6, PRKACB, and PRKCG (Figure 1D). In addition, to further verify the prognostic capacity of TRP channel-related genes, univariate Cox regression analysis was employed for screening (p < 0.05). TRPM7, MAPK11, and PRKCQ were identified as “protective” factors for ESCC patients, while PLA2G4D, RIPK3, CALML3 and MAPK9 were identified as “risk” factors (Figure 1E). Of note, in the TRP superfamily, only TRPM7 has a prognostic value in ESCC.

To explore new subtypes of TRP-related molecules in ESCC patients, we performed a consensus clustering analysis. In the GEO-ESCC cohort, the highest intra-group correlations and lowest inter-group correlations were observed when k = 2, indicating that patients could be divided into two clusters based on 18 TRP channel-associated DEGs (Figures 2A,B). Similarly, the TCGA-ESCC cohort was also divided into two clusters (Figures 2D,E). Survival analysis showed statistically significant differences between different molecular subgroups (p < 0.05, Figures 2C,F). There were statistically significant differences in staging and grading between subtypes in the GEO-ESCC cohort (Figure 2G). Staging between different subtypes in the TCGA-ESCC cohort also showed significant differences (Figure 2H). In summary, we have argued for the significance of TRP channel-related genes on survival and tumor progression in ESCC patients from another perspective.

To investigate the interaction between 18 TRP channel-associated DEGs, the online tool of the GeneMANIA website was employed to build the PPI network. Relatively close interactions were observed among the 18 TRP-channel-related DEGs (Figure 3A). We downloaded mutation data from TCGA database and analyzed 18 TRP channel-related DEGs. As shown in Figure 3B, missense mutations were the most common variant. Further, SNPs were the most common variant kind, and C > T was ranked as the top SNV class. Notably, TRPM3 had the highest mutation frequency among the 18 TRP channel-related DEGs.

We also revealed the CNV status of 18 TRP channel-related DEGs in ESCC patients. The results showed that most genes had heterozygous amplification and deletions (Figure 3C). Normal samples were used as controls to study the methylation of these genes. Interestingly, among the 18 TRP channel-related DEGs, only CALML3, PLA2G4E, PTGER2, HTR2C, and ADCY6 showed obvious diversities in the degree of methylation between normal samples and ESCC samples. Specifically, the degree of methylation of PTGER2, HTR2C, and ADCY6 was higher than that in normal samples; the degree of methylation of CALML3 and PLA2G4E was lower than that in normal samples (Figure 3D). Unfortunately, only the methylation degreels of above five TRP channel-related DEGs can be retrieved from the TCGA database. Finally, we estimated the role of each gene in tumor-related pathways using the GSVA algorithm, including apoptosis, DNA damage response, cell cycle, hormone AR, hormone ER, PI3K/AKT, EMT, RAS/MAPK, RTK, and TSC/mTOR. Based on the results, TRP channel-related genes can activate most of the above pathways, while inhibition only includes a few pathways, such as apoptosis and the cell cycle (Supplementary Figure S2).

Not only are ion channels dysregulated in cancer, but the expression of their regulators, effectors and other interacting genes is also significantly altered. TRP channels, as a typical ion channel, we speculate that the altered gene expression of related channels would reflect the risk of tumor progression from another perspective. To calculate TRP scores in ESCC patients, we screened 19 genes associated with prognosis in the GEO-ESCC cohort by univariate Cox regression analysis (p < 0.1), and further screened by LASSO regression analysis. The prognostic signature performed best when seven genes were included (Figures 4A,B). Finally, the regression coefficients of the seven genes were calculated using multivariate Cox regression analysis (Figure 4C). The signature formula is: TRP risk score = (−0.498711399 × expression level of TRPM6) + (−0.15894075 × expression level of HTR2C) + (−0.217314244 × expression level of PLA2G4A) + (−0.690911672 × expression level of ASIC4) + (−0.949344105 × expression level of MAPK14) + (−0.699670637 × expression level of PLCG2) + (0.673798617 × expression level of SRC).

Eighty patients with TGCA-ESCC were used as an external validation cohort. The gene expression data were standardised using the “sva” package before further analysis (Supplementary Figure S3). Based on the median value of risk scores, 80 patients in the TCGA cohort were divided into the high-risk group (n = 40) and the low-risk group (n = 40), and 119 patients in the GEO cohort were divided into the high-risk group (n = 59) and the low-risk group (n = 60). PCA analysis showed that both risk subgroups in both cohorts were well separated (Figures 5A,D). ROC analysis of the GEO-ESCC cohort showed a good predictive effect of the risk score (1-year AUC = 0.840, 2-year AUC = 0.785, and 3-year AUC = 0.794) (Figure 5B). Similarly, good predictive power was shown in the TGCA-ESCC cohort (1-year AUC = 0.813, 2-year AUC = 0.730, and 3-year AUC = 0.916) (Figure 5E). The Kaplan-Meier curve consistently showed that the OS of the high-risk subgroup was substantially shorter than that of the low-risk subgroup (p < 0.05, Figures 5C,F).

Combined with the clinical data, we again performed univariate and multivariate Cox regression analysis to assess whether the TRP risk score could be used as an independent predictor for ESCC. Univariate Cox regression analyses revealed that the risk score (p < 0.001, HR = 1.229), N stage (p = 0.002, HR = 1.443), and TNM stage (p = 0.004, HR = 1.901) were obviously associated with OS (Figure 6A). Multifactorial Cox regression analysis showed that only risk score was an independent predictor for OS (p < 0.001, HR = 1.200, Figure 6B).

Considering the importance of TNM staging in clinical practice, we combined TNM staging with TRP risk scores to construct a nomogram to predict 1-, 2-, and 3-year survival (Figure 6C). Besides, the calibration curve for the probability of 1, 2, and 3-year OS showed an optimal agreement between observation and prediction (Figure 6D). Moreover, TRP risk score with N stage (p = 0.044, Figure 6E), TNM stage (p = 0.017, Figure 6F), and T stage (p = 0.03, Figure 6G) showed significant differences.

These results suggest that TRP risk scores may influence the clinical outcomes of ESCC patients.

We further comprehensively evaluated the guiding role of TRP risk score in immunotherapy and immune-related mechanisms. As shown in Figures 7A,C, macrophages, such as those belonging to the M0, M1, and M2 subsets, constituted a significant proportion in both TCGA and GEO cohorts. The distribution of immune cell subsets between the low- and high-risk groups was shown using a box plot. In the low-risk group of the TCGA and GEO cohorts, the infiltrating scores of plasma cells was obviously higher and NK cells activated lower (Figures 7B,D). Meanwhile, in the GEO cohort, Pearson analysis was used to explore the relationship between seven signature-related genes and the level of 22 immune cell infiltrations. We found that ASIC4 was positively correlated with B cells memory, NK cells activated and negatively correlated with plasma cells. HTR2C was only positively correlated with T cells gamma delta. PLA2G4A was negatively correlated with macrophages M0, T cells CD8. Additional details were shown in Figure 7E. In addition, Figure 7F showed only the immune checkpoints with differential expression (TMIGD2, TNFRSF14, TNFRSF4, ICOS, CD80, LGALS9, HAVCR2, PDCD1, LAG3, LAIR1, CD40), all of which were upregulated in the high-risk group compared to the low-risk group.

The tumor purity score, estimate score, immune score, and stromal score were successfully generated with the ESTIMATE algorithm. Remarkably, patients who had a low-risk score demonstrated greater tumor purity (p = 0.038), lower estimate score (p = 0.038), lower immune score (p = 0.044), and higher stromal score (p = 0.045) than those who had a high-risk score (Supplementary Figures S4A–D), which was consistent with previous research findings that a lower estimate score indicates greater tumor purity. Human leukocyte antigen (HLA)-related genes exert a vital effect on controlling the immune response. After comparing the expression of HLA-associated genes among various groups, it was observed that most of the HLA-associated genes are upregulated in the high-risk group (Supplementary Figure S4E).

Additionally, the correlation between 7 signature-associated genes and immune infiltration in the ESCC microenvironment was established in order to get a better understanding of the relationship between these genes and immune infiltration. According to the findings, the expression of the majority of genes was shown to be negatively related to the level of immune cell infiltration (). Finally, we analyzed the transcriptional regulatory network of seven signature-associated genes, in which MAPK14 regulated the most TF and dominated the transcriptional regulatory network (). Supplementary Figure S5 Supplementary Figure S6

To further explore the potential mechanisms underlying the different clinical risks caused by abnormal TRP-related gene expression, we performed a GSEA analysis. The outcomes revealed that immune-associated pathways such as chemokine signaling pathway, cytokine-cytokine receptor interaction, adaptive immune response, and T cell activation were enhanced in the high-risk group (Figures 8A,B). Interestingly, the significant enrichment of metabolism- and cell growth-associated pathways was observed in the low-risk group, including arginine and proline metabolism, linoleic acid metabolism, terpenoid backbone biosynthesis, and keratinocyte differentiation (Figures 8C,D).

The current first-line chemotherapy regimens for ESCC in clinical guidelines include platinum-containing dual-agent chemotherapy regimens, paclitaxel combined with platinum regimens, or 5-fluorouracil-based regimens. In addition, anti-PD-1 drugs also have an important role in ESCC. We used the pRRophetic algorithm to evaluate the therapeutic effect of drugs. Unfortunately, there was no difference in PD-1 expression between different risk groups (Figure 9A), and correlation analysis also showed a weak correlation (Figure 9B). This reminded us that the novel immune checkpoint identified in Figure 7F might be useful for the treatment of patients with different TRP risks. In addition, the IC50 results were encouraging: the low-risk group showed better efficacy for paclitaxel (Figure 9C), cisplatin (Figure 9D), and 5-fluorouracil (Figure 9E) than the high-risk group, indicating that the TRP risk score plays a guiding role in the chemotherapy of ESCC patients.

We used the GTEx database and qRT-PCR experiment to analyze the expression levels of the 7 signature-associated genes were analyzed with the GTEx database and qRT-PCR experiment. As shown in Figures 10A–G, tumor tissues showed obviously higher expression levels of PLCG2, PLA2G4A, and HTR2C than the normal esophageal tissues. There was high expression of ASIC4 and MAPK14 in normal esophageal samples. Nevertheless, no obvious diversity was observed in the expression level of TRPM6 between ESCC and nearby tissues. The above-mentioned gene expression results in clinical tissue samples almost conformed to the RNA sequencing data analyzed by the TCGA database. In addition, results also showed that high-risk cohort had a shorter survival time in our hospital cohort (p < 0.05). Survival analysis also showed that high-risk cohort had a shorter survival time in our hospital cohort (p < 0.05) (Figures 10H).