Leveraging a disulfidptosis-based signature to improve the survival and drug sensitivity of bladder cancer patients

Background: Disulfidptosis is a recently discovered form of cell death. However, its biological mechanisms in bladder cancer (BCa) are yet to be understood.

Methods: Disulfidptosis-related clusters were identified by consensus clustering. A disulfidptosis-related gene (DRG) prognostic model was established and verified in various datasets. A series of experiments including qRT-PCR, immunoblotting, IHC, CCK-8, EdU, wound-healing, transwell, dual-luciferase reporter, and ChIP assays were used to study the biological functions.

Results: We identified two DRG clusters, which exhibited distinct clinicopathological features, prognosis, and tumor immune microenvironment (TIME) landscapes. A DRG prognostic model with ten features (DCBLD2, JAM3, CSPG4, SCEL, GOLGA8A, CNTN1, APLP1, PTPRR, POU5F1, CTSE) was established and verified in several external datasets in terms of prognosis and immunotherapy response prediction. BCa patients with high DRG scores may be characterized by declined survival, inflamed TIME, and elevated tumor mutation burden. Besides, the correlation between DRG score and immune checkpoint genes and chemoradiotherapy-related genes indicated the implication of the model in personalized therapy. Furthermore, random survival forest analysis was performed to select the top important features within the model: POU5F1 and CTSE. qRT-PCR, immunoblotting, and immunohistochemistry assays showed the enhanced expression of CTSE in BCa tumor tissues. A series of phenotypic assays revealed the oncogenetic roles of CTSE in BCa cells. Mechanically, POU5F1 can transactivate CTSE, promoting BCa cell proliferation and metastasis.

Conclusions: Our study highlighted the disulfidptosis in the regulation of tumor progression, sensitivity to therapy, and survival of BCa patients. POU5F1 and CTSE may be potential therapeutic targets for the clinical treatment of BCa.

1 Introduction

Bladder cancer (BCa) is the most common urological neoplasm, with a significant impact on public health worldwide. There were over 570,000 new cases and 210,000 deaths from BCa globally in 2020 (1, 2). BCa can be categorized into non-muscle-invasive and muscle-invasive BCa. While muscle-invasive BCa accounts for only about 30% of newly diagnosed cases, its aggressive nature, propensity for metastasis, drug resistance, and high rate of recurrence contribute to reduced cancer-specific survival after R0 resection (3, 4). Hence, elucidating the molecular mechanisms underlying the progression of BCa is of paramount importance.

Recently, Liu and colleagues uncovered a novel form of cell death, disulfidptosis, which has not been characterized previously (5). Disulfidptosis is induced by the accumulation of intracellular disulfides in glucose-starved cells with overexpressed SLC7A11 (6). Unlike ferroptosis and apoptosis, disulfidptosis is mediated by the susceptibility of the actin cytoskeleton to disulfide stress. The study also shows that glucose transporter inhibitors trigger disulfidptosis and control tumor proliferation, suggesting the significance of disulfidptosis in cancer management. Overall, the study sheds light on the role of disulfidptosis in controlling tumor progression.

In the study, we aimed to comprehensively investigate the role of DRGs in the prognosis, TIME landscapes, and drug resistance in BCa. We developed and validated a DRG-related prognostic model, which demonstrated high accuracy in predicting prognosis and response to immunotherapy across various independent cohorts. POU5F1 and CTSE were the top two important features within the model. Mechanistically, POU5F1 can transactivate CTSE and promote the proliferation and metastasis of BCa. These results suggested potential therapeutic targets for BCa treatment.

2 Methods

2.1 Data collection, tumor somatic mutation analysis, and protein-protein interaction analysis

The BCa datasets were obtained from the TCGA-BLCA and GEO databases (GSE13507, GSE32548, and GSE32894), as previously described (7). Additional immune checkpoint blockade (ICB) -treated datasets, including the Mariathasan BCa cohort (8), Braun RCC cohort (9), and Liu SKCM cohort (10), were obtained from original publications. The details of these included datasets were provided in Table S1. The details of DRGs were provided in Table S2.

Gene Ontology Biological Processes (GO-BP), KEGG, and cancer hallmarks were obtained from the MSigDB (https://www.gsea-msigdb.org/gsea/msigdb/) (11). Immunohistochemistry (IHC) images of proteins were procured from the Human Protein Atlas (https://www.proteinatlas.org/) (12).

The mutational landscape of all TCGA-BLCA patients and the mutation of DRGs were analyzed using the “maftool” package (13). Online tool GeneMANIA was employed to construct the PPI network of 14 DRGs (http://genemania.org/) (14).

2.2 Unsupervised consensus clustering, differential analysis, enrichment analysis, and TIME landscape estimation

Unsupervised consensus clustering analysis of the TCGA-BLCA cohort was conducted based on the expression profiles of 14 DRGs using the “ConsensusClusterPlus” package (15).

Differentially expressed genes (DEGs) between DRG-based molecular clusters were identified using the “limma” package (16). DEGs with P value (adjusted)< 0.01 and |logFC [fold change] | > 1 were considered significant. The details of DEGs were provided in Table S2.

The enrichment analysis (GO, KEGG, and gene set enrichment analysis [GSEA]) was employed via the “clusterProfiler” package as described previously (7, 17). The activity of GO-BP and KEGG terms for each sample was quantified by the “GSVA” package (18).

TIME scores including stromal, immune, and ESTIMATE scores were determined by the “ESTIMATE” package (19). The infiltration levels of 22 immune cell subsets were estimated by the “CIBERSORT” package (20).

2.3 Establishment and verification of a disulfidptosis−related prognostic model

The model was developed and validated as described in our previous study (21). In brief, the TCGA-BLCA dataset was considered as the training cohort, while the three GEO BCa datasets were external validation cohorts. Univariate Cox regression analysis was first used to identify overall survival (OS)-related DEGs (OS-DEGs). The details of OS-DEGs were provided in Table S2. The least absolute shrinkage and selection operator (LASSO) Cox regression analysis was then employed to reduce the dimensionality of high latitude data. The disulfidptosis-related predictive model was finally constructed using multivariate Cox regression analysis. Each BCa sample of TCGA-BLCA was assigned a DRG score using the formula: DRG score = ${\displaystyle \sum _{i=1}^{n}ki\ast Xi}$, where k_i represents the regression coefficient and X_i represents the relative expression level of gene i.

BCa samples were divided into high- and low-risk groups according to the median value of the DRG scores. Kaplan-Meier curves with the log-rank test were utilized to determine the OS and progression-free survival (PFS) between the two groups. The performance of the model was estimated by the time-dependent receiver operating characteristic (ROC) analysis.

R packages “survival”, “survminer”, and “timeROC” were used in these analyses (22, 23).

2.4 Tissue samples

The study protocol was approved by the Institutional Ethics Committee of Peking Union Medical College Hospital and we performed these experiments adhered to the principles of the Declaration of Helsinki. Tumor tissues of BCa and adjacent para-cancerous tissues were procured following radical surgery and histological confirmation. Written informed consent was obtained from all participants.

2.5 Cell culture and transfection

SV-HUC-1, T24, 5637, J82, and RT4 were purchased from the Cancer Institute of the Chinese Academy of Medical Sciences (Beijing, China). Cells were cultured in complete RPMI 1640 or DMEM supplemented with 10% fetal bovine serum (FBS) and 1% penicillin/streptomycin (Beyotime, Shanghai, China).

The siRNAs compounds targeting CTSE (siCTSE-1: 5’- CAACUACUUGGAUAUGGAAUA- 3’; siCTSE-2: 5’ – CAAUCUUUCUCCAUUCAGUAU – 3’) were designed by GenePharma (Suzhou, China). Transfection was performed using lipofectamine 3000 (Invitrogen) following the manufacturer’s instructions. Overexpression plasmid pcDNA3.1-POU5F1 and corresponding empty vector were purchased from Obio Technology Corp.

2.6 qRT-PCR, immunoblotting, and IHC

These assays were performed according to our previous studies (7, 17, 24). The primer sequences for qRT-PCR analysis were listed in Table S3. The information on all the antibodies used in the study was presented in Table S4.

2.7 Dual-luciferase reporter assay and chromatin immunoprecipitation

These procedures were performed as previously described (25). In brief, wild-type (CTSE-WT) and mutant-type (CTSE-Mut) sequences of the POU5F1 binding site in the promoter of CTSE were designed and cloned into the luciferase vectors (Beyotime, Shanghai, China). The overexpression vector or empty vector of POU5F1 was co-transfected into plated cells. After being cultured for 48 hours, dual luciferase activities were then evaluated by Dual-Luciferase Reporter Gene Assay Kit (Beyotime, Shanghai, China) following the manufacturer’s standard.

ChIP assay was employed using ChIP Assay Kit (Beyotime, Shanghai, China) following the instructions. After the cross-link, the anti-POU5F1 antibody was used to immunoprecipitate the corresponding genomic sequences. Then the sequences were analyzed by qPCR assay. The primer sequences for the CTSE promoter were provided in Table S1.

2.8 CCK-8, wound-healing, and transwell invasion assays

These phenotypic assays were performed as described in our previous studies (7, 17).

2.9 Statistical analyses

R 4.1.2 and GraphPad Prism 9.2.0 were employed for all statistical analyses. The correlation between the two variables was determined by the Spearman correlation coefficient. Student’s T-Test or Mann–Whitney U tests were used for continuous variables. Chi-Square Tests were utilized for categorized variables. A two-tailed p-value< 0.05 was considered statistically significant.

3 Results

3.1 Genetic and transcriptional characteristics of DRGs

First, we analyzed the somatic mutation landscapes of BCa patients, revealing that 94.69% of BCa samples showed somatic mutations (Figure 1A). Next, the somatic mutation pattern of 14 DRGs was also estimated (Figure 1B). Our results showed that DRGs in 112 out of 414 BCa samples (27.05%) had somatic mutations, primarily driven by missense mutations. In particular, MYH9, MYH10, ACTB, FLNB, and TLN1 showed the highest incidence of missense mutations, while MYH9 and ACTB had the highest incidence of nonsense mutations and in-frame deletions, respectively. To investigate the interactions between the 14 DRGs, the PPI network was constructed by the GeneMANIA online tool, which revealed that FLNA, MYH9, IQGAP1, CAPZB, and DSTN were regarded as the hub genes (Figure 1C).

FIGURE 1

Figure 1 Genetic and transcriptional characteristics of DRGs. (A) Somatic mutation landscape of all patients. (B) Somatic mutation of 14 DRGs. (C) PPI network of 14 DRGs by GeneMANIA. (D) mRNA levels of 14 DRGs. ns, not significance, *p< 0.05, ***p< 0.001, ****p< 0.0001. ns, not significance.

We also analyzed the expression profiles of the 14 DRGs in the TCGA-BLCA cohort. Results showed that ACTB, DSTN, FLNA, IQGAP1, MYL6, and TLN1 were overexpressed in normal tissues, while CD2AP and INF2 were overexpressed in tumor tissues (Figure 1D). We further found that ACTB, DSTN, FLNA, and TLN1 were upregulated in high-grade and high-stage (stage III and IV) BCa tissues, indicating their association with BCa progression (Figures S1A, B). Finally, representative IHC images demonstrated the protein levels of the four genes (Figure S1C).

3.2 DRG-based molecular clusters with distinct clinical features and TIME landscapes

To investigate the biological roles of DRGs in BCa, we employed an unsupervised consensus clustering algorithm to categorize BCa samples of the TCGA-BLCA cohort based on the expression profiles of 14 DRGs. As indicated by the CDF curves and the PAC test, the optimal clustering number was 2 (Figures S2A–C). The expression levels of DRGs were higher in BCa samples of C1 compared to those of C2 (Figure 2A). In terms of clinical features, the proportion of patients with age > 65, male gender, and stage III and IV was relatively higher in C1 than in C2 (Figures 2A; S2D–F). Moreover, patients in C1 had a poorer prognosis compared to those in C2 (Figure 2B).

FIGURE 2

Figure 2 DRG-based molecular clusters with distinct clinical features and TIME landscapes. (A) Expression profiles of DRGs and clinicopathological characteristics between clusters. (B) Survival analysis between C1 and C2. (C) Expression levels of CRGs between clusters. (D) Expression profiles of ICGs and clinicopathological characteristics between clusters. (E) Differences in TIME scores between clusters. (F) Abundances of infiltrating immune cells between clusters. ns, not significance, *p< 0.05, **p< 0.01, ***p< 0.001, ****p< 0.0001. ns, not significance.

Next, we evaluated the expression levels of immune checkpoint genes (ICGs, Table S5) and chemoradiotherapy-related genes (CRGs, Table S5), as well as the TIME scores between two DRG-based clusters. As presented in Figure 2C, SOX2, EZH2, EGFR, and AKR1C1 of the CRGs had distinct expression patterns between the two DRG-based clusters. Specifically, AKR1C1 and SOX2 were downregulated in C1, whereas EGFR and EZH2 were overexpressed in C1. Regarding the TIME landscapes, we observed that samples of C1 had overexpressed ICGs (Figure 2D), a higher TIME score (Figure 2E), and higher infiltration levels of macrophages and activated memory CD4+ T cells (Figure 2F). However, samples of C2 had higher abundances of regulatory T cells, activated dendritic cells, monocytes, infiltrated memory B cells, and neutrophils.

3.3 DRGs-based molecular clusters with dysregulated pathways and biological process

To investigate the pathways and biological processes between the two DRG clusters, a series of enrichment analyses were employed, including GSVA, GSEA, and over-representation analysis (ORA). The results of GO-BP GSVA showed that DRG C1 was enriched in cellular structure-related processes, including podosome assembly, membrane raft assembly, regulation of protein maturation, and cortical cytoskeleton organization (Figure 3A). KEGG term results revealed that DRG C1 was abundant in immunity-related pathways, including focal adhesion and NOD-like receptor signaling pathway (Figure 3B). GSEA subsequently uncovered that DRG C1 was significantly linked to cancer hallmarks, including the cell cycle (G2M checkpoint and mitotic spindle, Figure 3C) and cancergenesis and progression (epithelial-mesenchymal transition [EMT], hypoxia, and angiogenesis, Figure 3D). In contrast, DRG C2 was associated with metabolism, specifically oxidative phosphorylation and fatty acid metabolism (Figure 3E).

FIGURE 3

Figure 3 DRGs-based molecular clusters with dysregulated pathways and biological process. (A, B) GSVA of GO-BP (A) and KEGG (B) terms between clusters. (C–E) GSEA of significant hallmarks enriched in DRG C1 (C, D) and C2 (E). (F, G) GO (F) and KEGG (G) analysis of cluster DEGs.

We then identified DRGs cluster-related DEGs to verify the findings above. Consistently, cellular structure and immunity-related biological processes and pathways were mainly enriched (Figures 3F, G).

3.4 Detection of gene clusters related to disulfidptosis in BCa

First, prognosis-related DEGs were screened out by univariate Cox regression analysis. A hierarchical clustering algorithm was then used to classify the BCa samples into three gene clusters (Figure 4A). We observed that DRG C1 was highly correlated with gene cluster 1 while C2 was mostly related to gene cluster 3. Consistently, BCa patients in gene cluster 1 had a worse prognosis compared to those in gene cluster 3 (Figure 4B). Additionally, the majority of DRGs had significantly varied expression levels among the three gene clusters (Figure 4C).

FIGURE 4

Figure 4 Detection of gene clusters related to disulfidptosis. (A) The expression profiles of DEGs and the distribution of DRG clusters among gene clusters 1 to 3. (B) Survival analysis among three gene cluster. (C) mRNA levels of 14 DRGs among gene clusters 1 to 3. (D) Distribution of DRG scores between DRG clusters. (E) mRNA levels of 14 DRGs between groups. (F) Distribution of DRG scores among three gene clusters. ns, not significance, **p< 0.01, ***p< 0.001, ****p< 0.0001.

3.5 Development and validation of a disulfidptosis-related prognostic model

Based on the OS-related DEGs, a disulfidptosis-related prognostic model was established for potential clinical application. The TCGA-BLCA cohort was used as the training set whereas three GEO datasets (GSE13507, GSE32548, and GSE32894) were regarded as the testing sets. Followed by LASSO and multivariate Cox regression analysis, a model with ten genes was constructed, including DCBLD2, JAM3, CSPG4, SCEL, GOLGA8A, CNTN1, APLP1, PTPRR, POU5F1, CTSE (Figure S3). The DRG score of each BCa sample was: DRG score =0.03458897* DCBLD2 + 0.03469218* JAM3 + 0.05195818* CSPG4 + 0.04814716* SCEL - 0.14422347* GOLGA8A + 0.02423207* CNTN1 + 0.09220036* APLP1 - 0.05143026* PTPRR + 0.02817304* POU5F1 + 0.03568603* CTSE.

Next, we aimed to investigate the correlation between DRG clusters (C1 and C2), gene clusters (1, 2, and 3), and DRG scores. Our findings indicated that DRG scores were considerably higher in DRG C1 than in C2. Moreover, the expression of DRGs was higher in BCa patients from the high-risk group (Figures 4D, E). We also discovered that the DRG score in gene clusters followed a rank order of 1 > 2 > 3, as illustrated in Figure 4F.

As illustrated in the risk plot, the DRG score was linked to higher mortality and decreased survival (Figure 5A). BCa patients of the high-risk group had poorer OS than those of the low-risk group (Figure 5B). Time-dependent ROC analysis showed excellent performance of the model in predicting the prognosis: AUC of 0.837, 0.833, and 0.83 at 1-, 3-, and 5-year time points, respectively (Figure 5C).

FIGURE 5

Figure 5 Development and validation of the disulfidptosis-related prognostic model. (A) Upper panel: Survival time and status between risk groups. The red dots represent dead BCa patients whereas the others indicate the alive. Bottom panel: Distribution of the DRG scores in TCGA-BLCA cohort. (B) Survival analysis between risk groups. (C) Prognostic performance of the model. (D, F, H) Survival analysis between risk groups in the GSE13507 (D), GSE32548 (F), and GSE32894 (H). (E, G, I) Prognostic performance of the model in GSE13507 (E), GSE32548 (G), and GSE32894 (I).

To further verify the performance, we tested it on three external validation cohorts. Similarly, a declined survival was observed in BCa patients of the high-risk group (Figures 5D, F, H). High AUC value from time-dependent ROC analysis also suggested the prognostic performance of the established model (Figures 5E, G, I).

3.6 Estimation of DRG score in immune infiltration

The biological features of BCa samples in high- and low-risk groups were assessed by GSEA (Figures S4A, B). The analysis revealed that cellular structure-related activities such as external encapsulating structure organization and collagen fibril organization, as well as inflammation-related GO-BP including leukocyte migration, inflammatory response, and neutrophil chemotaxis, were highly activated in the high-risk group. Conversely, metabolism-related biological processes, such as long-chain fatty acid and arachidonic acid metabolic processes, were suppressed in the high-risk group (Figure S4A). Furthermore, the results of GSEA on cancer hallmarks indicated that a high DRG score was linked to cell cycle (G2M checkpoint and E2F target), tumor progression (EMT, hypoxia, and angiogenesis), as well as inflammation and immunity (inflammatory response and TNFA signaling via NF-kappa B) (Figure S4B).

TIME scores including immune, stromal, and ESTIMATE scores were remarkably elevated in the high-risk group (Figure S4C). The correlation between immune infiltration and ten model genes was assessed (Figure S4D). Macrophages showed a positive correlation with DRG score, whereas CD8+ T cells, follicular helper T cells (Tfh cells), and Tregs exhibited a negative correlation with DRG score (Figure S4E).

Cancer stem cells possess self-renewal, multipotent, and tumor-initiating abilities, leading to tumor growth, recurrence, and resistance to current treatments. To evaluate the relationship between the DRG score and the CSC index, both were combined, revealing a weak negative correlation (Figure S4F).

3.7 Correlations between DRG score and ICGs and tumor mutational burden

Patients in the high-risk group (52%) exhibited a higher frequency of TP53 mutations than those in the low-risk group (45%) (Figures S5A, B). Additionally, Patients in the high-risk group had significantly higher TMB scores compared to those in the low-risk group (Figure S5C).

Given the reported association between ICB therapy and the expressions of ICGs, we next evaluated the association between the ICGs and DRG score. Our findings indicated a significant negative correlation between most ICGs and the ten model genes (Figure S5D). Specifically, the expressions of PVR, CD276, and SIRPA increased as the DRG score increased, while the expressions of TNFRSF14 and CD96 had the opposite trends (Figure S5E).

3.8 Evaluation of the model in ICB-treated cohorts

Our previous findings demonstrated a potential correlation between DRG scores and the expression of ICGs. To further investigate the predictive ability of our developed model in immunotherapy response, we analyzed three ICB-treated cohorts comprising different types of tumors. Within the Mariathasan cohort, the high-risk group was associated with a worse prognosis (Figure 6A). Non-responders had significantly higher risk scores compared to responders (Figure 6B), and the high-risk group had a lower proportion of responders and a higher proportion of non-responders than the low-risk group (Figure 6C). The ROC analysis revealed a high predictive performance of the model (Figure 6D). Additionally, the risk score gradually decreased from the desert immune phenotype to the inflamed phenotype (Figure 6E).

FIGURE 6

Figure 6 The performance of the model in predicting immunotherapeutic response in ICB-treated cohorts. (A) Survival analysis between risk groups. (B) Distribution of the DRG score between responders and non-responders. (C) Distribution of the responders and non-responders between risk groups. (D) The immunotherapeutic response prediction performance of the model. (E) Distribution of the DRG score among three immune phenotypes. (A–E) Data was analyzed in the Mariathasan cohort. (F, K, G, I) The Kaplan–Meier OS (F, K) and PFS (G, I) curves between risk groups. (H, M) Distribution of the DRG score between responders and non-responders. (I, N) Distribution of the responders and non-responders between high- and low-risk group. (J, O) The immunotherapeutic response prediction performance of the model. (F-J) Data was analyzed in the Braun RCC cohort. (K-O) Data was analyzed in the Liu SKCM cohort.

Within the Braun cohort, renal cell carcinoma (RCC) patients in the high-risk group suffered from worse OS and PFS (Figures 6F, G). Non-responders had significantly higher risk scores compared to responders (Figure 6H), and the high-risk group had a lower proportion of responders and a higher proportion of non-responders compared to the low-risk group (Figure 6I). The ROC analysis uncovered a high predictive performance of the model (Figure 6J). Similar findings were observed in the Liu cohort of patients with skin cutaneous melanoma (SKCM) (Figures 6K–O).

3.9 Correlation between DRG score and CRGs and chemotherapeutic sensitivity

Except for MGMT and AKR1C1, a negative correlation was observed between the expressions of CRGs and the ten model genes (Figure S6A). Additionally, EZH2, EGFR, TBX5, and SOX2 were positively correlated with the DRG score (Figure S6B). These findings suggested that the DRG score may be able to predict the chemotherapeutic response. To verify the findings, drug sensitivity analysis was performed using six frequently used chemotherapeutic agents in BCa patients. The results demonstrated that BCa patients in the high-risk group exhibited resistance to cisplatin, doxorubicin, gemcitabine, methotrexate, paclitaxel, and vinblastine (Figure S6C).

3.10 POU5F1 transactivates CTSE by directly binding to its promoter

The random survival forests (RSF) algorithm was utilized to ascertain the relative importance of each gene in the model. As presented in Figure 7A, POU5F1 displayed the highest degree of significance, followed by CTSE, DCBLD2, and GOLGA8A. Our correlation analysis indicated that POU5F1 was positively correlated exclusively with CTSE (R = 0.39, Figure 7B). Furthermore, this result was consistent with the outcomes obtained through analysis of the GEPIA online tool (Figure 7C). POU5F1, also known as OCT4, is a transcriptional factor that has been linked with tumor proliferation, migration, and therapy resistance (26). Previous research has thoroughly explored the oncogenic phenotypes of POU5F1 in various human cancers, including BCa (27). Additionally, Fristrup et al. examined the protein levels of CTSE (Cathepsin E) in a large, multicenter cohort and found that CTSE levels were significantly associated with progression to stage T2-T4 BCa (28). Nonetheless, the functional roles of CTSE in BCa remain poorly understood. In light of the positive correlation between POU5F1 and CTSE, we hypothesized that POU5F1 might transactivate CTSE to promote BCa progression.

FIGURE 7

Figure 7 POU5F1 promoted BCa proliferation and metastasis by transactivating CTSE. (A) The relative importance of ten model genes. (B) Correlation among the ten model genes. (C) Spearman correlation between POU5F1 and CTSE by GEPIA. (D, E) Relative expression levels of CTSE in BCa tumor and adjacent normal tissues. (F). Representative IHC images of CTSE in BCa tissues and adjacent normal tissues. (G) Expression levels of CTSE in several BCa cell lines by qRT-PCR and immunoblotting. (H) Immunoblotting and qRT-PCR validated the knockdown efficacy of siRNAs targeting CTSE in 5637 cell line. (I, J) Deficient CTSE inhibited 5637 proliferation as indicated by CCK-8 (I) and EdU (J) assays. (K, L) Deficient CTSE inhibited 5637 migration (K) and invasion (L). (M) Diagram illustrated the predicted binding site of POU5F1 to CTSE promoter. (N) Luciferase activity of the POU5F1/CTSE promoter reporter was examined. (O) Luciferase assay with different doses of POU5F1 overexpression plasmid. (P) ChIP analysis indicated the enrichment of POU5F1 on the gene promoter region of CTSE. (Q) qRT-PCR and immunoblotting detected CTSE expression following POU5F1 depletion. (R) Rescuing CCK-8 assays. (S) Rescuing wound-healing assays. (T) Rescuing transwell invasion assays. (N-T) Assays were performed in 5637 cells. *p< 0.05, **p< 0.01, ***p< 0.001.

First, we investigated the expression levels of CTSE in clinical BCa tissues and BCa cell lines. In 106 pairs of BCa tissues, we observed an upregulation of both CTSE mRNA and protein levels in tumor tissues (Figures 7D, E), which was further supported by representative IHC images (Figure 7F). In particular, CTSE was found to be overexpressed in BCa cell lines, especially in 5637 (Figure 7G). To gain insights into the biological functions of CTSE in BCa, a series of phenotypic assays were conducted in 5637 cells. We verified the efficacy of CTSE knockdown using qRT-PCR and immunoblotting assays (Figure 7H). CTSE deficiency was found to inhibit 5637 cell line proliferation (Figures 7I, J), as well as wound-healing migration and transwell invasion assays (Figures 7K, L). Overall, the upregulation of CTSE drives the progression of BCa cells in vitro.

Next, we investigated the regulatory mechanisms between the transcription factor POU5F1 and the gene CTSE. A previous study reported the existence of a transcriptional binding site for POU5F1 (29). We then identified a potential binding site for POU5F1 at the promoter of CTSE (Figure 7M). To investigate the functional significance of this binding site, we conducted luciferase assays in 5637 cells. We found that co-transfection of POU5F1 significantly stimulated the luciferase activity of the CTSE-WT promoter (Figure 7N). Moreover, increasing the dose of POU5F1 resulted in a gradual enhancement of luciferase activity (Figure 7O). Next, a ChIP assay was employed to validate the enrichment of POU5F1 on the promoter of CTSE (Figure 7M) and found that exogenous POU5F1 promoted more POU5F1 enrichment on the CTSE promoter compared to the vector group (Figure 7P). Our results demonstrate that POU5F1 transactivates CTSE by directly binding to its promoter, which has regulatory effects on both mRNA and protein levels (Figure 7Q).

3.11 POU5F1 promotes BCa proliferation and metastasis by transactivating CTSE

We examined the functional association between POU5F1 and CTSE through phenotypic assays as POU5F1 was capable of regulating CTSE by physically binding to its promoter. The CCK-8 assay demonstrated that overexpression of POU5F1 enhanced the proliferation of 5637, while deficiency of CTSE inhibited it (Figure 7R). Similar outcomes were observed in wound-healing migration (Figure 7S) and transwell invasion assays (Figure 7T). In summary, our findings suggested that POU5F1 boosts the proliferation and metastasis of BCa cells by transactivating CTSE.

4 Discussion

Cell death is a crucial factor in regulating tumor proliferation (30). Previous research has established a strong association between cell death and cancer cell metabolism (31). However, the underlying mechanisms linking cell death and metabolism in BCa remain poorly understood.

Disulfidptosis, initially proposed by Xiaoguang Liu, has provided new insights into the role of disulfides and glucose metabolism dysregulation in cell death. However, the landscapes of DRGs in BCa remain unclear. In this study, we conducted a systematic investigation of the genetic and transcriptional changes of 14 DRGs in BCa. We also established a disulfidptosis-related prognostic model with 10 features which exhibited excellent performance in predicting prognosis and immunotherapeutic response.

By deciphering the genetic and transcriptional landscapes of DRGs in BCa, we found that only 27.05% of BCa samples had genetic mutations in DRGs. Surprisingly, six DRGs (ACTB, DSTN, FLNA, IQGAP1, MYL6, and TLN1) were downregulated in BCa compared to normal tissues, while CD2AP and INF2 were upregulated. Most of these downregulated DRGs (except for IQGAP1) were overexpressed in high-grade tumor tissues. Additionally, mRNA levels of ACTB, DSTN, FLNA, and TLN1 increased gradually with the tumor stage. Besides, protein levels of these genes were higher in BCa tumor tissues than in normal tissues, indicating post-translational modifications (32–35).

Based on the expression pattern of DRGs, two disulfidptosis-related molecular clusters of the TCGA-BLCA cohort were identified. BCa patients in DRG C1 were characterized by decreased survival and advanced clinicopathological features. Further decoding of the TIME unraveled that C1 was featured by the inflamed TIME in terms of upregulated expression profiles of ICGs, elevated TIME scores, and infiltrated immune cell subsets. As evidenced by previous BCa studies, tumors with distinct TIME landscapes may hold different sensitivities to chemotherapy and immunotherapy (7, 17, 36, 37). Consistently, the expression levels of CRGs and ICGs varied significantly between DRG clusters, indicating varied therapeutic responses between DRG clusters. These findings indicated the importance of disulfidptosis in driving BCa.

Further, a disulfidptosis-related prognostic model was developed based on the OS-DEGs of the DRG clusters. Deciphering the TIME of BCa unveiled that a high DRG score was linked to the inflamed phenotype which had significant effects on the cancer treatment. Furthermore, a higher RNA stemness score was detected in the high-risk group, indicating the crucial role of disulfidptosis patterns in maintaining BCa tumors.

Using the RFS method, we identified an exclusively positive correlation between two critical genes, POU5F1 and CTSE. POU5F1, also known as OCT4, encodes a transcription factor regulating stem cell pluripotency via the POU homeodomain. POU5F1 regulates the characteristics of tumor-initiating cells in terms of survival, self-renewal, resistance to drugs, and EMT (38). Roundhill et al. studied the role of POU5F1 in Ewing sarcoma and found that POU5F1 intensely interacted with stemness and chemoresistance genes (39). Mitchell and colleagues addressed the critical roles of POU5F1 in WDR5-induced glioblastoma progression (40). As for CTSE (Cathepsin E), its biological characteristics have been addressed in various tumors including pancreatic cancer (41), prostate cancer (42), and gastric cancer (43). In our study, we observed higher expression levels of CTSE at both mRNA and protein levels in BCa tumor tissues. Moreover, overexpression of CTSE promoted the proliferation, migration, and invasion of BCa cells. Mechanistically, we found that POU5F1 directly binds to the promoter of CTSE, leading to its transactivation and promotion of BCa progression.

5 Conclusion

In the study, we comprehensively investigated the DRG profiles in BCa and established a disulfidptosis-related prognostic model which exhibited excellent performance in predicting prognosis and immunotherapeutic response. BCa sample of different DRG scores was also characterized by distinct TIME landscape, response to immunotherapy and chemotherapy, and dysregulated pathways and biological processes. Furthermore, we found that POU5F1 and CTSE were critical components of the prognostic model. We also conducted further investigations to uncover the regulatory mechanisms underlying the relationship between POU5F1 and CTSE.

Data availability statement

This data can be found here: TCGA-BLCA: https://portal.gdc.cancer.gov/GSE13507: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE13507 GSE32548: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE32548 GSE32894: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE32894.

Ethics statement

The studies involving human participants were reviewed and approved by Institutional Ethics Committee of Peking Union Medical College Hospital. The patients/participants provided their written informed consent to participate in this study.

Author contributions

Conceptualization, HC and ZJ. Methodology, HC. Software, HC. Validation, HC. Formal analysis, HC. Investigation, HC. Resources, WY. YL, and LM. Data curation, WY, YL, and LM. Writing—original draft preparation, HC. Writing—review and editing, HC. Visualization, HC. Supervision, WY, YL, LM, and ZJ. Project administration, ZJ. All authors contributed to the article and approved the submitted version.

Funding

This study was supported by the National High Level Hospital Clinical Research Funding (2022-PUMCH-B-008).

Acknowledgments

We thank Dr. Jianming Zeng (University of Macau), and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes.

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.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2023.1198878/full#supplementary-material

Supplementary Figure 1 | Transcriptional and protein changes of DRGs. (A) mRNA levels of DRGs between low- and high-grade BCa tumor tissues. (B) mRNA levels of DRGs among stage II, III, and IV. (C) The protein levels of ACTB, DSTN, FLNA, and TLN1 between tumor and normal tissues by the Human Protein Atlas.

Supplementary Figure 2 | DRG-based molecular clusters with distinct clinicopathological features. (A) The consensus score matrix when clustering number was 2. (B) The CDF curves of each clustering number. (C) The proportion of ambiguous clustering (PAC) score indicated the optimal clustering number of 2. (D–F) The distribution of age groups (> 65 and ≤ 65) (D), gender (E), and stage (F) between DRG-based clusters.

Supplementary Figure 3 | LASSO regression.

Supplementary Figure 4 | The correlation between the DRG score and TIME score. (A, B) The activated or suppressed GO-BP (A) and hallmarks (B) in DRG C1. (C) The correlations between TIME scores and DRG scores. (D) The correlations between the ten model genes and infiltration of 22 immune cell subsets. (E, F) Correlations between the DRG score and the abundances of immune cells (E) and the stem index (F).

Supplementary Figure 5 | Genetic variations and ICGs. (A, B) The mutation features of BCa patients in low- (A) and high-risk (B) groups. (C) Distribution of TMB between risk groups. (D) Correlations between ICGs and ten model genes. (E) Correlations between ICGs and the DRG score.

Supplementary Figure 6 | The correlation between DRG score and CRGs. (A) The correlations between ten model genes and CRGs. (B) The correlations between CRGs and the DRG score. (C) Predicted IC50 between risk groups.

Supplementary Table 1 | The details of included datasets.

Supplementary Table 2 | The details of the 14 DRGs, DEGs between DRG-based molecular clusters, and OS-DEGs.

Supplementary Table 5 | The lists of ICGs and CRGs.

References

1. Chen H, Yang W, Xue X, Li Y, Jin Z, Ji Z. Neoadjuvant immunotherapy and chemoimmunotherapy for stage II-III muscle invasive bladder cancer. Front Immunol (2022) 13:986359. doi: 10.3389/fimmu.2022.986359

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Lenis AT, Lec PM, Chamie K, Mshs MD. Bladder cancer: a review. Jama. (2020) 324(19):1980–91. doi: 10.1001/jama.2020.17598

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Malmström PU, Rintala E, Wahlqvist R, Hellström P, Hellsten S, Hannisdal E. Five-year followup of a prospective trial of radical cystectomy and neoadjuvant chemotherapy: Nordic cystectomy trial i. the Nordic cooperative bladder cancer study group. J Urol. (1996) 155(6):1903–6. doi: 10.1016/S0022-5347(01)66042-7

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Chen X, Zhang J, Ruan W, Huang M, Wang C, Wang H, et al. Urine DNA methylation assay enables early detection and recurrence monitoring for bladder cancer. J Clin Invest. (2020) 130(12):6278–89. doi: 10.1172/JCI139597

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Liu X, Nie L, Zhang Y, Yan Y, Wang C, Colic M, et al. Actin cytoskeleton vulnerability to disulfide stress mediates disulfidptosis. Nat Cell Biol (2023) 25(3):404–414. doi: 10.1038/s41556-023-01091-2

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Koppula P, Zhuang L, Gan B. Cystine transporter SLC7A11/xCT in cancer: ferroptosis, nutrient dependency, and cancer therapy. Protein Cell (2021) 12(8):599–620. doi: 10.1007/s13238-020-00789-5

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Ma L, Chen H, Yang W, Ji Z. Crosstalk between mesenchymal stem cells and cancer stem cells reveals a novel stemness-related signature to predict prognosis and immunotherapy responses for bladder cancer patients. Int J Mol Sci (2023) 24(5):4760. doi: 10.3390/ijms24054760

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Mariathasan S, Turley SJ, Nickles D, Castiglioni A, Yuen K, Wang Y, et al. TGFβ attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. Nature. (2018) 554(7693):544–8. doi: 10.1038/nature25501

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Braun DA, Hou Y, Bakouny Z, Ficial M, Sant’ Angelo M, Forman J, et al. Interplay of somatic alterations and immune infiltration modulates response to PD-1 blockade in advanced clear cell renal cell carcinoma. Nat Med (2020) 26(6):909–18. doi: 10.1038/s41591-020-0839-y

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Liu D, Schilling B, Liu D, Sucker A, Livingstone E, Jerby-Arnon L, et al. Integrative molecular and clinical modeling of clinical outcomes to PD1 blockade in patients with metastatic melanoma. Nat Med (2019) 25(12):1916–27. doi: 10.1038/s41591-019-0654-5

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Liberzon A, Subramanian A, Pinchback R, Thorvaldsdóttir H, Tamayo P, Mesirov JP. Molecular signatures database (MSigDB) 3.0. Bioinformatics (2011) 27(12):1739–40. doi: 10.1093/bioinformatics/btr260

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Uhlén M, Fagerberg L, Hallström BM, Lindskog C, Oksvold P, Mardinoglu A, et al. Proteomics. tissue-based map of the human proteome. Science. (2015) 347(6220):1260419. doi: 10.1126/science.1260419

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Mayakonda A, Lin DC, Assenov Y, Plass C, Koeffler HP. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res (2018) 28(11):1747–56. doi: 10.1101/gr.239244.118

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Warde-Farley D, Donaldson SL, Comes O, Zuberi K, Badrawi R, Chao P, et al. The GeneMANIA prediction server: biological network integration for gene prioritization and predicting gene function. Nucleic Acids Res (2010) 38:W214–20. doi: 10.1093/nar/gkq537

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics. (2010) 26(12):1572–3. doi: 10.1093/bioinformatics/btq170

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res (2015) 43(7):e47–e. doi: 10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Chen H, Yang W, Xue X, Li Y, Jin Z, Ji Z. Integrated analysis revealed an inflammatory cancer-associated fibroblast-based subtypes with promising implications in predicting the prognosis and immunotherapeutic response of bladder cancer patients. Int J Mol Sci (2022) 23(24):15970. doi: 10.3390/ijms232415970

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinf (2013) 14(1):7. doi: 10.1186/1471-2105-14-7

CrossRef Full Text | Google Scholar

19. Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun (2013) 4:2612. doi: 10.1038/ncomms3612

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods (2015) 12(5):453–7. doi: 10.1038/nmeth.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Chen H, Pan Y, Jin X, Chen G. An immune cell infiltration-related gene signature predicts prognosis for bladder cancer. Sci Rep (2021) 11(1):16679. doi: 10.1038/s41598-021-96373-w

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Blanche P, Dartigues J-F, Jacqmin-Gadda H. Estimating and comparing time-dependent areas under receiver operating characteristic curves for censored event times with competing risks. Stat Med (2013) 32(30):5381–97. doi: 10.1002/sim.5958

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Lin H, Zelterman D. Modeling Survival Data: Extending the Cox Model. Technometrics (2002) 44(1):85–86. doi: 10.1198/tech.2002.s656

CrossRef Full Text | Google Scholar

24. Yang W, Chen H, Ma L, Dong J, Wei M, Xue X, et al. A comprehensive analysis of the FOX family for predicting kidney renal clear cell carcinoma prognosis and the oncogenic role of FOXG1. Aging (2022) 14(24):10107–24. doi: 10.18632/aging.204448

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Liu W, Shen D, Ju L, Zhang R, Du W, Jin W, et al. MYBL2 promotes proliferation and metastasis of bladder cancer through transactivation of CDCA3. Oncogene. (2022) 41(41):4606–17. doi: 10.1038/s41388-022-02456-x

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Hallis SP, Kim SK, Lee JH, Kwak MK. Association of NRF2 with HIF-2α-induced cancer stem cell phenotypes in chronic hypoxic condition. Redox Biol (2023) 60:102632. doi: 10.1016/j.redox.2023.102632

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Mao M, Yang L, Hu J, Liu B, Liu C, Zhang X, et al. OCT3/4 enhances tumor immune response by upregulating the TET1-dependent NRF2/MDM2 axis in bladder cancer. Genomics. (2021) 113(4):2122–33. doi: 10.1016/j.ygeno.2021.04.033

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Fristrup N, Ulhøi BP, Birkenkamp-Demtröder K, Mansilla F, Sanchez-Carbayo M, Segersten U, et al. Maspin, Plk1, and survivin are promising prognostic protein markers for progression in non-muscle invasive bladder cancer. Am J Pathol (2012) 180(5):1824–34. doi: 10.1016/j.ajpath.2012.01.023

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Rodda DJ, Chew J-L, Lim L-H, Loh Y-H, Wang B, Ng H-H, et al. Transcriptional regulation of nanog by OCT4 and SOX2*. J Biol Chem (2005) 280(26):24731–7. doi: 10.1074/jbc.M502573200

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Bao J-h, Lu W-c, Duan H, Ye Y-q, Li J-b, Liao W-t, et al. Identification of a novel cuproptosis-related gene signature and integrative analyses in patients with lower-grade gliomas. Front Immunol (2022) 13. doi: 10.3389/fimmu.2022.933973

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Villa E, Ricci JE. How does metabolism affect cell death in cancer? FEBS J (2016) 283(14):2653–60. doi: 10.1111/febs.13570

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Varland S, Vandekerckhove J, Drazic A. Actin post-translational modifications: the Cinderella of cytoskeletal control. Trends Biochem Sci (2019) 44(6):502–16. doi: 10.1016/j.tibs.2018.11.010

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Treppiedi D, Catalano R, Mangili F, Mantovani G, Peverelli E. Role of filamin a in the pathogenesis of neuroendocrine tumors and adrenal cancer. Endocrine Oncol (2022) 2(1):R143–R52. doi: 10.1530/EO-22-0055

CrossRef Full Text | Google Scholar

34. Malla RR, Vempati RK. Talin: a potential drug target for cancer therapy. Curr Drug Metab (2020) 21(1):25–32. doi: 10.2174/1389200221666200214114018

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Tan L, Cheng W, Liu F, Wang DO, Wu L, Cao N, et al. Positive natural selection of N6-methyladenosine on the RNAs of processed pseudogenes. Genome Biol (2021) 22(1):180. doi: 10.1186/s13059-021-02402-2

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Zhang Q, Liu S, Wang H, Xiao K, Lu J, Chen S, et al. ETV4 mediated tumor-associated neutrophil infiltration facilitates lymphangiogenesis and lymphatic metastasis of bladder cancer. Adv Sci (2023) 10(11):e2205613. doi: 10.1002/advs.202205613

CrossRef Full Text | Google Scholar

37. Xie R, Cheng L, Huang M, Huang L, Chen Z, Zhang Q, et al. NAT10 drives cisplatin chemoresistance by enhancing ac4C-associated DNA repair in bladder cancer. Cancer Res (2023) 83(10):1666–83. doi: 10.1158/0008-5472.22746976

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Wang Y-J, Herlyn M. The emerging roles of Oct4 in tumor-initiating cells. Am J Physiology-Cell Physiol (2015) 309(11):C709–C18. doi: 10.1152/ajpcell.00212.2015

CrossRef Full Text | Google Scholar

39. Roundhill EA, Pantziarka P, Liddle DE, Shaw LA, Albadrani G, Burchill SA. Exploiting the stemness and chemoresistance transcriptome of Ewing sarcoma to identify candidate therapeutic targets and drug-repurposing candidates. Cancers (2023) 15(3):769. doi: 10.3390/cancers15030769

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Mitchell K, Sprowls SA, Arora S, Shakya S, Silver DJ, Goins CM, et al. WDR5 represents a therapeutically exploitable target for cancer stem cells in glioblastoma. Genes Dev (2023) 37(3-4):86–102. doi: 10.1101/gad.349803.122

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Ye H, Li T, Wang H, Wu J, Yi C, Shi J, et al. TSPAN1, TMPRSS4, SDR16C5, and CTSE as novel panel for pancreatic cancer: a bioinformatics analysis and experiments validation. Front Immunol (2021) 12:649551. doi: 10.3389/fimmu.2021.649551

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Yasukochi A, Kawakubo T, Nakamura S, Yamamoto K. Cathepsin e enhances anticancer activity of doxorubicin on human prostate cancer cells showing resistance to TRAIL-mediated apoptosis. Biol Chem (2010) 391(8):947–58. doi: 10.1515/bc.2010.087

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Konno-Shimizu M, Yamamichi N, Inada K, Kageyama-Yahara N, Shiogama K, Takahashi Y, et al. Cathepsin e is a marker of gastric differentiation and signet-ring cell carcinoma of stomach: a novel suggestion on gastric tumorigenesis. PloS One (2013) 8(2):e56766. doi: 10.1371/journal.pone.0056766

PubMed Abstract | CrossRef Full Text | Google Scholar