Using Necroptosis-Associated Genes To Predict The Immune Microenvironment And Prognosis Of Bladder Urothelial Carcinoma

doi:10.21203/rs.3.rs-2073511/v1

Download PDF

Research Article

Using Necroptosis-Associated Genes To Predict The Immune Microenvironment And Prognosis Of Bladder Urothelial Carcinoma

https://doi.org/10.21203/rs.3.rs-2073511/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

PURPOSE

Bladder urothelial carcinoma (BLCA), the most common urinary tract malignancy, has a high recurrence rate and poor survival at late stages. Necroptosis, a form of programmed cell death, is involved in cancer development and progression, but its function in BLCA prognosis remains unclear. This study sought to investigate the role of necroptosis in the development and prognosis of BLCA.

METHODS

Clinical information and RNA expression matrix data were obtained from the databases. Survival analysis was performed to obtain survival- and necroptosis-related genes and identify any that overlapped. Consensus clustering analysis was used to create different subgroups by combining the overlapping gene expression matrix and clinical information. The tumor immune microenvironment and immune status of the different subgroups were determined using ESTIMATE, MCPcounter, and ssGSEA analysis. We performed differential analysis on the gene expression matrix of molecular subpopulations to find and screen out differentially expressed genes (DEGs). GO, KEGG, GSVA, and GSEA analyses were used to elucidate the underlying mechanisms of the DEGs. Lasso Cox regression analysis was used to build a prognostic risk model and perform a pan-cancer analysis of the screened genes. The results were used to define potential roles for these genes in other cancers and assess the efficacy of the risk model.

RESULTS

Cluster analysis identified two subgroups, C1 and C2, with significantly different survival rates. ESTIMATE, MCPcounter, and ssGSEA analyses showed that high immune scores, tumor purity, and immune status were associated with poorer prognoses. GO and KEGG functional enrichment analyses indicated that DEGs were mainly focused on tumor proliferation, invasion, and immunity and GSEA analysis suggested that necroptosis may affect Toll-like receptor signaling pathways, MAPK cascade regulation of leukocyte trafficking, and cytokine-cytokine receptor interaction pathways. Lasso Cox regression analysis was used to model the prognostic risk while screening for representative necroptosis-associated genes, ANXA1, ATAD3A, and TRPC6, with high potential for survival prediction in BLCA patients. The pan-cancer analysis indicated that the three representative genes were also differentially expressed in other cancer types.

CONCLUSION

Expression of necroptosis-related genes such as ANXA1, ATAD3A, and TRPC6 correlate with the immune microenvironment of BLCA patients and have the potential for use in disease prognostics.

necroptosis

bladder urothelial carcinoma

immune microenvironment

prognosis

individualized therapy

pan-cancer analysis

Bladder urothelial carcinoma (BLCA) is the fourth most common cancer in men and the most commonly diagnosed urological malignancy worldwide(Green et al., 2013). The treatment of BLCA is closely related to the stage of tumor progression(Witjes et al., 2014;Patel et al., 2020). Once BLCA has developed a muscular infiltrate, not only is the 5-year survival rate less than 40%, but there is an increased risk of recurrence or progression, and even after radical surgery, 50% of patients still have recurrence after surgery(Witjes et al., 2014;Patel et al., 2020). Most recurrences occur when distant metastases have already occurred, resulting in metastatic bladder cancer. The current standard of care for muscle-invasive bladder cancer (MIBC) is radical cystectomy followed by neoadjuvant chemotherapy, but the quality of life for patients after surgery is poor(Smith et al., 2013). In addition, approximately 30% of patients with BLCA are diagnosed with muscle-invasive bladder cancer (MIBC) or metastatic disease(Giridhar and Kohli, 2017). To improve the prognosis of advanced BLCA, there is a critical need for individualized and targeted therapeutic strategies, risk stratification methods, and the identification of new prognosis-related genes.

Necroptosis is a newly identified form of programmed cell death(Gong et al., 2019). Several studies suggest that necroptosis is involved in cancer development and tumor responsiveness to treatment(Gong et al., 2019). For example, findings indicate that inducing necroptosis can help to control breast cancer progression(Thakur et al., 2019). In addition, combination berberine and cisplatin treatment induces necroptosis and apoptosis in ovarian cancer cells(Liu et al., 2019). Thus, targeting necroptosis-related pathway genes is being investigated as a new therapeutic strategy for cancer treatment(Della Torre et al., 2021). However, the role of necroptosis in BLCA and its ability to predict disease outcomes remains unclear.

This study comprehensively analyzed necroptosis- and prognosis-related genes in BLCA to explore how necroptosis impacts the tumor immune microenvironment (TIME) and the survival of BLCA patients. A risk score model was constructed to assess the prognostic value of necroptosis-related genes in BLCA, providing new insight into the underlying molecular mechanisms of this disease and the development of targeted therapeutic strategies. Meanwhile, We performed a pan-cancer analysis of the screened genes to give some hints for their prognosis in the development of other cancers.

2.1 Patients and Datasets

Clinical, prognostic, and sequenced RNA-seq data from standardized BLCA patient samples were downloaded from The Cancer Genome Atlas website (TCGA) (https://portal.gdc.cancer.gov/) and defined as the training cohort. BLCA patient data was also obtained from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/; GSE13507) for use as the validation cohort. Information on genes associated with necroptosis was obtained from the GENECARDs (https://www.genecards.org/) database.

2.2 Screening for prognosis-related genes

The prognostic significance of each gene was assessed using the Cox Proportional-Hazards model in the R “survivor” package(Simon et al., 2011). Survival time, survival status, and gene expression data were integrated from clinical samples (P < 0.05)(Tibshirani, 1997).

2.3 Cluster analysis

Genes obtained from the survival analysis were compared with genes associated with necroptosis to identify any that overlapped. The overlapping genes were clustered using ConsensusClusterPlus(Wilkerson and Hayes, 2010;Yu et al., 2012). 80% of the clusters were sampled 10 times in duplicate. An empirical cumulative distribution function plot was used to determine the optimal number of clusters(Rooney et al., 2015) .

2.4 TIME evaluation and immune infiltration analysis

The number of mesenchymal and immune cells in malignant tissues was estimated using the ESTIMATE to calculate the mesenchymal score, immune score, and tumor purity(Li et al., 2020). The MCP counter was used to calculate the abundance of ten immune infiltrating cells, CD4 + T cells, CD8 + T cells, cytotoxic lymphocytes, B cells, NK cells, monocytes, myeloid dendritic cells, neutrophils, endothelial cells, and fibroblasts by the R “ggExtra,” “vioplot,” “ggpubr,” and “gglot2” package(Yoshihara et al., 2013) .Single sample gene set enrichment analysis (ssGSEA) was performed using the R package by calculating the enrichment fraction of each sample and gene set pair and determining the enrichment of each immune cell type.

2.5 Functional enrichment analysis

The R “limma” package was used to screen for differentially expressed genes (DEGs) between the two molecular subgroups. The DEGs were subjected to Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses and related pathways were identified using the R “ClusterProfiler” package and visualized in Metascape (https://metascape.org/gp/index.html)(Yu et al., 2012). The enrichment score was determined for each sample in the expression matrix using Gene Set Variation Analysis (GSVA) (DOI:10.18129/B9.bioc.GSVA, version 1.40.1) and used to rank each of the genes. Relevant pathways and molecular mechanisms were predicted from the gene expression profiles(Walter et al., 2015).

2.6 Risk model building and validation

The R “glmnet” package was used to integrate survival time, survival status, and gene expression data for regression analysis using the Lasso Cox regression method. In particular, a 10-fold cross-validation was used to obtain the optimal model for the genes. The risk score for each patient in the training and validation sets was calculated as follows:

RiskScore = 0.185996363365641 × TRPC6 expression value + 0.107698904611713 × ANXA1 expression value + 0.208717048548769 × ATAD3A expression value.

The predictive efficiency of the model was assessed using the receiver operating characteristic (ROC) method.

2.7 Pan-cancer analysis

A high-quality prognostic dataset was obtained from previously published prognostic studies of TCGA using the UCSC Cancer Browser (https://xenabrowser.net/datapages/)(Chen et al., 2021). TARGET follow-up data were obtained as a supplement for samples with a follow-up time of < 30 days. The log2 (x + 0.001) transformation was further applied to each expression value and cancers with < 10 samples in a single cancer species were excluded to obtain final overall survival data for multiple cancer species. Pan-cancer prognosis data for the optimal model genes were analyzed using the R package.

2.8 Data analysis

The data were all statistically analyzed using R (version 3.6.0) and GraphPad Prism (version 8.0). Survival analysis was performed using the Kaplan–Meier method and the predictive performance of the risk model was assessed using the time-dependent ROC of the R “SurvivalROC” package. Subgroup analysis was performed by age, sex, BMI, tumor status, and clinical stage. Discontinuous data were expressed as numbers and percentages, while continuous data were expressed as the mean ± standard deviation (SD). Statistical analysis between two groups was performed using the t-test, with a flexible choice of one-way ANOVA when three or more groups were present. Differences were statistically significant when p < 0.05.

3.1 Clustering by necroptosis-related genes identified two molecular subtypes

Clinical information and RNA expression profile matrices for 412 samples were used as a training cohort and survival analysis was performed on the samples (Table 1), and the whole process of data analysis was depicted (Fig. 1). Results identified 7,171 genes associated with BLCA patient prognosis (Log rank test; P-value < 0.05) (Supplementary Datasheet S1). A total of 614 genes associated with necroptosis were obtained from GENECARDs (Supplementary Datasheet S2) and 108 overlapping genes associated with BLCA prognosis and necroptosis were screened (Fig. 2A, Supplementary Datasheet S3). Cluster analysis was performed using ConsensusClusterPlus and clusters were evaluated based on their average intra-cluster concordance. The number of clusters with the highest average intra-group consistency and best aggregation stability was 2 (k = 2,185 and 227 for C1 and C2, respectively) (Fig. 2B–D, Supplementary Figure S1-S2). Overlapping gene expression in the two subtypes was determined from a heat map and significant differences were found for most overlapping genes (Fig. 2E). To better resolve the subtypes, Kaplan-Meier curve analysis was performed using clinical information, and C1 patients had a higher overall survival than C2 patients (P < 0.001) (Fig. 2F). These findings verified that overlapping BLCA prognosis- and necroptosis-related genes classified BLCA patients into one of two molecular subtypes: 1) C1, with higher overall survival and 2) C2, with lower overall survival.

Table 1

**The clinical features in the training and verification cohort.**
Characteristics	Training(N = 412)	Verification(N = 165)	Total(N = 577)	pvalue	FDR
Age
Mean ± SD	68.10 ± 10.58	65.18 ± 11.97	67.26 ± 11.07
Median[min-max]	69.00[34.00,90.00]	66.00[24.00,88.00]	68.00[24.00,90.00]
Sex				0.05	0.05
Female	108(18.72%)	30(5.20%)	138(23.92%)
Male	304(52.69%)	135(23.40%)	439(76.08%)
BMI
Mean ± SD	27.09 ± 6.23	Na	27.09 ± 6.23
Median[min-max]	26.00[14.96,68.32]	Na	26.00[14.96,68.32]
Stage				1.80E-121	9.10E-121
Stage I	2(0.35%)	Na	2(0.35%)
Stage II	131(22.70%)	Na	131(22.70%)
Stage III	141(24.44%)	Na	141(24.44%)
Stage IV	136(23.57%)	Na	136(23.57%)
Unknown	2(0.35%)	165(28.60%)	167(28.94%)
Grade				3.80E-52	1.10E-51
High Grade	388(67.24%)	60(10.40%)	448(77.64%)
Low Grade	21(3.64%)	105(18.20%)	126(21.84%)
Unknown	3(0.52%)	0	3(0.52%)
Recurrence				1.50E-38	2.90E-38
No	246(42.63%)	67(11.61%)	313(54.25%)
Unknown	0	62(10.75%)	62(10.75%)
Yes	166(28.77%)	36(6.24%)	202(35.01%)
Cancer_Status				8.80E-92	3.50E-91
TUMOR FREE	235(40.73%)	0	235(40.73%)
Unknown	40(6.93%)	165(28.60%)	205(35.53%)
WITH TUMOR	137(23.74%)	0	137(23.74%)
Invasiveness				5.10E-126	3.60E-125
Unknown	412(71.40%)	0	412(71.40%)
muscle invasive	0	62(10.75%)	62(10.75%)
non-muscle invasive	0	103(17.85%)	103(17.85%)
Systemic chemo				5.10E-126	3.60E-125
No	0	138(23.92%)	138(23.92%)
Unknown	412(71.40%)	0	412(71.40%)
Yes	0	27(4.68%)	27(4.68%)

3.2 Immunological analysis

The tumor immune microenvironment (TIME) is created by the close relationship between tumorigenesis, growth, and metastasis and the internal and external environment in which tumor cells are located(Pitt et al., 2016). Tumor-infiltrating immune cells are the primary non-tumor component of this environment and play an important role in predicting cancer patient prognosis(Yoshihara et al., 2013). Thus, ESTIMATE, MCPcounter, TIMER and ssGSEA analyses were performed on patients with both subtypes. The immune status of each group was assessed using the R “Estimate” package, and the immune score (P < 0.0001), assessment score (P < 0.0001), and tumor purity (P < 0.0001) were significantly higher in patients with BLCA subtype C2 than those with C1 (Fig. 3A). The MCP counter algorithm showed that B cells (P < 0.0001), cytotoxic lymphocytes (P < 0.0001), monocytes (P < 0.0001), myeloid dendritic cells (P < 0.0001), endothelial cells (P < 0.001), fibroblasts (P < 0.0001), NK cells (P < 0.0001), dendritic cells (P < 0.0001), and T cells (P < 0.0001) were significantly more abundant in patients with the C2 than C1 subtype (P < 0.0001) (Fig. 3B). The TIMER also suggests similar results (Supplementary Figure S3). In ssGSEA, all immune cell types were significantly higher in patients with the C2 subtype than those with the C1 subtype (Fig. 3C), and ssGSEA algorithm-generated immune cell lineage signature genes were significantly differentially expressed between patients with the C1 and C2 subtypes, with a relatively low expression in those with C1 (Fig. 3D). These results suggest that there are significant differences in the TIME and immune status between patients with these two subtypes; those with C2 had a higher ESTIMATE score, increased enrichment of genes associated with the immune cell lineage, and higher numbers of immune cells than those with C1.

3.3 DEG expression and functional analyses

Differential analysis of the C1 and C2 expression matrices was conducted using the R “limma” package to finally obtain DEGs. Based on a difference multiplicity of 2, P-value < 0.05, and FDR < 0.05, 1,021 DEGs were screened, of which 593 were downregulated in patients with the C2 subtype and 428 were upregulated in those with the C1 subtype (Fig. 4A, B, Supplementary Datasheet S4). The R “Profiler” package was used to perform enrichment analysis of the DEGs, including GO, KEGG, GSVA, and GSEA. GO enrichment analysis indicated the DEGs participate in many molecular functions including identical protein binding, signaling receptor binding, structural molecule activity, calcium ion binding, matrix structural constituent, glycosaminoglycan binding, integrin binding, growth factor binding, and collagen binding (Fig. 4C). The DEGs were mainly enriched in biological processes such as locomotion, biological adhesion, cell migration, regulation of immune system process, response to biotic stimulus, cell activation, response to cytokine, regulation of cell adhesion, external encapsulating structure organization (Fig. 4D). The DEGs were mainly involved in cellular components such as the external encapsulating structure, cell surface, collagen containing extracellular matrix, secretory vesicle, secretory granule, anchoring junction, side of membrane, endoplasmic reticulum lumen, collagen trimer (Fig. 4E). KEGG enrichment analysis identified several pathways associated with tumor proliferation, invasion, and immunity, including cell adhesion molecule camouflage, focal adhesion, ECM receptor interactions, bladder cancer development, cytokine-cytokine receptor interactions, antigen processing and presentation, Toll-like receptor signaling, Nod-like receptor signaling, and pathways involved in cancer (Fig. 5A).

To further explore the relationship between the enriched pathways and BLCA patient prognosis, GSVA and GSEA analyses were performed to assess differences in the relative expression of pathways regulated by the two clusters. GSVA analysis revealed that patients with the C2 subtype had higher enrichment in the Fc gamma R-mediated phagocytosis, cell adhesion molecules cams, Nod-like receptor signaling pathway, JAK/STAT signaling pathway, cytokine-cytokine receptor interactions, chemokine signaling pathway, Toll-like receptor signaling pathway, focal adhesion, ECM receptor interaction than those with the C1 subtype (Fig. 5B). GSEA analysis showed that Toll-like receptor signaling pathway, MAPK cascade regulating leukocyte trafficking, and cytokine-cytokine receptor interactions were highly expressed in patients with the C2 subtype (Fig. 5C–E). These results suggest that genes with differential C1 and C2 expression are closely associated with tumor proliferation, differentiation, invasion, and immune dysfunction in BLCA patients and may also be linked to poor disease prognosis.

3.4 Training queue risk modeling

The training cohort was used to extract overlapping necroptosis- and prognosis-related genes and construct a risk model for assessing BLCA patient prognosis. Lasso Cox analysis was used to evaluate the overlapping gene expression matrix and 44 potential genes with the best lambda values were screened (Table 2, Fig. 6A, B). Multivariate Cox analysis was performed using qualifying hazard ratios (HR) > 1 and P-values < 0.05, and three genes, Annexin-A1 (ANXA1), ATPase Family AAA Domain Containing 3A (ATAD3A), and transient receptor potential cation channel protein 6 (TRPC6), were identified and used to construct the risk model (Fig. 6C, D). The model successfully classified the BLCA patients into high- and low-risk groups and Multivariate Survival Analysis indicated that these three genes may be associated with BLCA patient prognosis. Notably, all three candidate genes were more highly expressed in the high- than the low-risk group and the overall survival rate of patients in the high-risk group was poorer than those in the low-risk group (Fig. 6E–G). ROC analysis was performed to verify the accuracy of the model. The constructed risk model had accurate predictive power at 5 years, with an area under the curve for the 1-, 3- and 5-year ROC curves of 0.80, 0.83, and 0.84, respectively (Fig. 6H). ESTIMATE analysis was also performed to verify the efficacy of the high- and low-risk subgroups for the TIME assessment. The stromal Score (p < 0.001) and ESTIMATE Score (p = 0.02) were significantly higher in the high- than the low-risk group, while immune purity were not statistically significant (p > 0.05) (Fig. 7A). ESTIMATE results also demonstrated that the higher immune score in the high-risk group was associated with poorer disease prognosis. Thus, the constructed risk model was correlated significantly with TIME and reliably predicted BLCA patient outcomes.

Table 2

**44 potential genes by Lasso Cox analysis.**
Gene	Coefficients	HR	95%CI Lower	95%CI Upper	pvalue
HNRNPH3	-0.60384694	0.359374814	0.228150391	0.566075108	1.01E-05
H1-2	-0.085387265	0.852833698	0.759110021	0.958128989	0.00736066
SP6	-0.121905942	0.836976347	0.73215614	0.956803293	0.00913946
ANXA1	0.107698905	1.177532673	1.029304148	1.347107362	0.017278206
ATAD3A	0.208717049	1.637993327	1.076741092	2.49179878	0.021144724
GNLY	-0.15791446	0.798835843	0.655147166	0.974038715	0.02642173
TRPC6	0.185996363	1.694493519	1.036912309	2.769094611	0.035324812
LACTB	-0.30258027	0.631230406	0.396498983	1.004925214	0.052470275
ITPK1	-0.154766098	0.717754656	0.510296079	1.009554586	0.056737222
MAPK3	0.136083693	1.368599692	0.978166909	1.914872708	0.067079165
TYRO3	0.128415358	1.359018307	0.950445275	1.943226829	0.092687214
DIRAS3	0.135061772	1.315032715	0.946500809	1.82705712	0.102625861
G6PD	0.055718668	1.149397251	0.96729209	1.365786047	0.113626042
PADI4	0.078433084	1.329967856	0.923178925	1.916003984	0.125806689
AURKC	-0.131139386	0.760385241	0.518492691	1.115127995	0.160861979
EGR1	0.023978805	1.093095201	0.951271481	1.256063219	0.209330788
VDAC2	0.115829173	1.251745975	0.873128799	1.79454393	0.221801567
CALM1	0.186276705	1.254799445	0.837748357	1.879468498	0.270847434
SIRT6	-0.153189454	0.777733387	0.474671532	1.274290072	0.318371786
CDK9	-0.075161125	0.796113558	0.507430696	1.249031251	0.321068117
IKBKB	-0.185846632	0.830045769	0.570997634	1.206617923	0.329099295
AIM2	-0.036090507	0.941152152	0.828595081	1.068999072	0.350684191
CFLAR	-0.050589204	0.812532578	0.507159814	1.301777414	0.387986493
GSK3A	-0.053571871	0.839735225	0.535875485	1.315893836	0.445973067
MLKL	-0.063391118	0.870802089	0.589885568	1.285497255	0.486332157
TNIP2	-0.020601632	0.866359325	0.567930127	1.321603563	0.505539799
FNDC4	0.003304879	1.061107583	0.887393609	1.268827375	0.51552998
LRRC59	0.08075975	1.162140141	0.736165529	1.834600583	0.518888741
CCT6A	0.023959957	1.140375215	0.763051592	1.704282705	0.521667885
LRP1	0.038298552	1.054247599	0.820810822	1.354073279	0.679109273
KRT1	0.006994068	1.017684573	0.935356623	1.10725884	0.683792642
CLTC	0.089164953	1.105946401	0.664777015	1.839891291	0.698194527
FASN	0.081399829	1.040522025	0.826109962	1.310583498	0.735816324
TRAF5	-0.083039055	0.934946953	0.601458976	1.453342356	0.76504448
YTHDC1	-0.10963455	0.93419186	0.526352787	1.658040869	0.816103814
ASIC1	0.127394623	1.037488829	0.717857106	1.49943918	0.844719578
LTBP1	0.009375392	0.980893228	0.808402304	1.190188993	0.844999389
GJA1	0.000753684	0.990344789	0.871850145	1.124944244	0.881380813
PPP1R12C	-0.019979834	0.972994363	0.639209302	1.481076741	0.898376543
EGFR	0.013720382	0.9888699	0.831503752	1.17601836	0.899285726
MAP1B	0.070334363	1.015307169	0.799566251	1.28925983	0.900806707
IGF2BP3	0.007071624	0.989961628	0.787542678	1.244407512	0.931112359
AKNA	-8.46039E-05	0.986425354	0.688844492	1.412561167	0.940530157
CASP6	-0.018380697	0.988391179	0.716703039	1.363070992	0.943235903

3.5 Screening for independent predictors of BLCA patient prognosis

To identify independent predictors of BLCA patient prognosis, the relationship between risk scores and patient clinical characteristics (age, sex, clinical stage, tumor status, BMI and tumor malignancy) was explored and risk models were constructed and evaluated for each characteristic. The risk scores of each sample were obtained from the Lasso Cox analysis in relation to their corresponding clinical characteristics. Age > 60 years, male gender, tumor status, and more advanced cancer stage were associated with higher risk scores (P < 0.05), while obesity (BMI > 30) had no significant effect on risk scores (Fig. 7B–G). Thus, the age, gender, and tumor status variables were used to regroup BLCA patients as follows: male BLCA group, female BLCA group, BLCA group > 60 years of age, BLCA group ≤60 years of age, with tumor BLCA group, and tumor-free BLCA group. KM curve analysis and risk model construction were conducted for each group. The risk model still had strong predictive power and a lower risk score remained associated with a better prognosis (Fig. 8A–F). Moreover, univariate/multivariate Cox regression analysis showed that the constructed risk model was an independent predictor of BLCA patient prognosis (Tables 3,4).

Table 3

**Univariate cox analysis in the training and verification cohort.**
Variates	HR	(95%CI)		Likelihood ratio test		Wald test		logrank test
Variates	HR	Lower	Upper	W value	p value	W value	p value	W value	p value
Training
Tumor Status	0.24	0.17	0.34	76.18	2.6E-18	72.61	1.6E-17	84.87	3.2E-20
Stage	0.58	0.48	0.71	31.51	0.00000002	29.87	4.6E-08	31.25	2.3E-08
Tumor Malignancy	2.86	0.71	11.57	3.14	0.08	2.18	0.14	2.39	0.12
BMI	0.78	0.53	1.13	1.81	0.18	1.73	0.19	1.74	0.19
age	2.08	1.43	3.02	16.76	0.000042	14.6	0.00013	15.26	0.000094
Sex	1.13	0.82	1.57	0.55	0.46	0.56	0.45	0.56	0.45
Riskscore	5.17	3.74	7.15	107.67	3.2E-25	99.09	2.4E-23	119.82	6.9E-28
Verification
Recurrence	7.7E-10	0	NaN	6.21	0.01	0	1	5.44	0.02
Systemic chemo	4.25	2.1	8.63	13.77	0.00021	16.08	0.000061	19.06	0.000013
Sex	2.1	0.97	4.54	3.14	0.08	3.54	0.06	3.7	0.05
Tumor Malignancy	5.99	2.76	12.99	23.58	0.0000012	20.47	0.000006	26.32	2.9E-07
age	4.92	2.35	10.27	19.19	0.000012	17.95	0.000023	21.87	0.0000029
Invasiveness	24.74	7.51	81.54	54.78	1.3E-13	27.81	1.3E-07	59.8	1E-14
Riskscore	2.52	1.23	5.15	5.91	0.02	6.4	0.01	6.86	0.0088

Table 4

**Multivariate cox analysis in training and verification cohort.**
Multivariate Cox analysis	Coefficient	exp(coef)	se(coef)	z	exp(coef)	exp(-coef)	lower 0.95	upper0.95	p-value	Concordance	Likelihood ratio test	Wald test	logrank test
Training										0.797 (se = 0.017 )	p = 2E-16	p = 2E-16	p = 2E-16
sex	0.007043	1.007068	0.180979	0.039	1.0071	0.993	0.7063	1.4358	0.969
age	0.021862	1.022103	0.008706	2.511	1.0221	0.9784	1.0048	1.0397	0.012
Clinical	0.131424	1.140452	0.224509	0.585	1.1405	0.8768	0.7345	1.7708	0.558
Tumor_status	-1.110115	0.329521	0.266834	-4.16	0.3295	3.0347	0.1953	0.5559	3.18E-05
grade	0.580185	1.786368	0.743885	0.78	1.7864	0.5598	0.4157	7.6766	0.435
recurrence	-0.193351	0.824193	0.255243	-0.758	0.8242	1.2133	0.4998	1.3592	0.449
riskscore	1.24977	3.489539	0.135643	9.214	3.4895	0.2866	2.6749	4.5523	2.00E-16
Verification										0.902 (se = 0.016)	p = 2E-14	p = 2E-14	p = 2E-14
sex	8.27E-01	2.29E + 00	4.33E-01	1.909	2.29E + 00	4.37E-01	0.9781	5.348	0.05628
age	5.27E-02	1.05E + 00	1.92E-02	2.742	1.05E + 00	9.49E-01	1.0151	1.095	0.00611
invasiveness	2.12E + 01	1.55E + 09	5.55E + 03	0.004	1.55E + 09	6.48E-10	0	Inf	0.99696
grade	2.24E-01	1.25E + 00	4.34E-01	0.515	1.25E + 00	8.00E-01	0.5343	2.926	0.60646
recurrence	1.91E + 01	2.05E + 08	5.55E + 03	0.003	2.05E + 08	4.89E-09	0	Inf	0.99725
riskscore	6.24E-01	1.87E + 00	2.99E-01	2.091	1.87E + 00	5.36E-01	1.04	3.351	0.0365

3.6 BLCA risk model relationship to the TIME and disease prognosis

The clinical information and RNA expression matrix of the BLCA samples obtained from the GEO database were used as a validation cohort to further evaluate the established prognostic risk model and determine its relationship to the TIME. BLCA patients in the validation cohort were classified as high or low risk using the formula described above to derive corresponding risk scores. Expression of the three genes, ANXA1, ATAD3A, and TRPC6, previously screened by heat map was higher in the high-risk group (Fig. 9A). Survival analysis of patients in the high- and low-risk groups indicated that those in the high-risk group had a poorer prognosis (p < 0.001) (Fig. 9B). ROC analysis showed that the risk model had the highest predictive accuracy in predicting 5-year survival (Fig. 9C). The relationship between the risk model and the TIME was also assessed. The results matched what was observed for the training cohort, with significantly higher interstitial, immune, and estimation scores in the high-risk than in the low-risk group, while tumor purity remained equivalent (Fig. 9D). These results suggest that the developed risk model was associated with both disease prognosis and the TIME in the BLCA validation cohort.

3.7 Construction and calibration of a BLCA integrated model

An integrated model that combined both risk models and clinical characteristics was constructed to more accurately predict BLCA patient prognosis. Survival time, survival status, and seven characteristics, age, sex, tumor stage, tumor grade, recurrence, tumor status, and risk score were integrated to build a nomogram using Cox regression and assign it a specific score (Fig. 10A). The nomogram was validated in both the training and validation cohorts (Supplementary Figure S4). The nomogram, C-index, and calibration curve were used to diagnose the composite model and demonstrated acceptable accuracy. The overall C-index of the model was: 0.7995, 95% CI (0.7649–0.8341). The calibration curves showed that the actual survival rates at 1,3 and 5 years were well matched in the validation and training cohorts (Fig. 10B–G). The results indicated that the integral nomogram could accurately predict BLCA patient prognosis.

BLCA is the most common malignant tumor of the urinary system and one of the 10 most common cancers(Bray et al., 2018). The high recurrence rate and low prognosis of advanced BLCA make treatment a challenge. This study analyzed and predicted BLCA patient prognosis and related target genes based on their clinical characteristics. Findings should inform the development of more individualized and targeted treatment options for this disease.

Several recent studies have identified new methods for predicting BLCA patient prognosis. For example, miR-133 downregulation was shown to predict the progression and poor prognosis of patients with urothelial carcinoma of the bladder(Zhao et al., 2013), and let-7c, miR-100, and miR-145 expression were associated with the EMT score and abundance of macrophages in tumor tissue(Feng et al., 2019). In addition, stromal lymphocytic infiltrate (SLI) was shown to be a reliable marker of poor BLCA patient prognosis and positively correlated with histological grade, tumor stage, and lymph node status(Furuya et al., 2019). The presence of CD163 + and CD103 + macrophages also predicted the poor prognosis of patients with primary T1 high-grade bladder urothelial carcinoma(Yang et al., 2019;Jin et al., 2022). However, the predictors identified to date remain insufficient for reliable BLCA prognosis and precision treatment.

To determine more reliable prognosis markers, survival analysis using clinical information and the RNA matrix of BLCA patients was conducted to screen for prognosis-related genes and compare them with programmed necroptosis-related genes to identify genes that overlapped. Different molecular subtypes were obtained using consensus clustering analysis of the combined overlapping gene expression matrix along with clinical information. Immunological analysis of the molecular subgroups was also performed to determine the TIME and immune status of the different subtypes. Furthermore, differential analysis of gene expression matrices of the molecular subtypes was performed to find and screen out DEGs and functional analysis was used to identify the underlying mechanisms. Finally, Lasso and multifactorial Cox regression analyses were used to build a prognostic risk model and screen out representative genes. A pan-cancer analysis of the screened genes was then used to characterize the potential mechanisms of their gene targets on other cancer prognoses and to evaluate the efficacy of the risk model.

Lasso regression analysis was used to screen three genes, TRPC6, ATAD3A, and ANXA1, for prognostic relevance based on their role in necroptosis, which was shown to be closely associated with tumor development. TRPC6, a member of the TRPC family, can affect breast cancer proliferation and migration through Ca2⁺-dependent signaling, while inhibition of TRPC6 reduces proliferation and invasion of non-small cell lung cancer cells(Yang et al., 2017). LncRNA TUG1 promotes cell growth and migration in colorectal cancer through the TUG1-miR-145-5p-TRPC6 pathway(Wang et al., 2021). Helicobacter pylori can promote gastric cancer migration and invasion by upregulating TRPC6 through Wnt/β-catenin signaling(Song et al., 2019).

The AAA structural domain of ATAD3A, a nuclear-encoded mitochondrial enzyme, is involved in cell death processes. Transcriptional upregulation and increased protein stability are potential mechanisms for elevated ATAD3A in cancer(Teng et al., 2019). ATAD3A also mediates the activation of RAS-independent mitochondrial ERK1/2 signaling, which facilitates the development of head and neck cancer. ANXA1 acts as a damage-associated molecular pattern-linked molecule that promotes necroptosis(Lang et al., 2022). By binding to and stabilizing EphA2, ANXA1 induces the growth and metastasis of nasopharyngeal carcinoma. Transcriptome analysis was used to confirm ANXA1 as a prognostic and immune microenvironmental marker for glioma(Feng et al., 2020). ANXA1-containing extracellular vesicles also regulate macrophage polarization in the TIME and promote the development and metastasis of pancreatic cancer. And the typical immunostaining graphs of three genes in BLCA tumor tissue were provided by the human protein atlas (HPA) database (Supplementary Figure S5)(Novizio et al., 2021).

The current study used pan-cancer analysis to determine the prognosis-related expression of these three representative genes in other tumors. Specifically, the R “survival” package was used to build a risk regression model and analyze the prognostic relationship between gene expression and specific cancers by using the Log rank test. It was observed that TRPC6 prognosed poor outcomes when highly expressed in tumor models such as glioma (TCGA-GBMLGG), stomach and esophageal carcinoma (TCGA-STES), kidney renal papillary cell carcinoma (TCGA-KIRP), stomach adenocarcinoma (TCGA-STAD), bladder urothelial carcinoma (TCGA-BLCA), mesothelioma (TCGA-MESO) and acute lymphoblastic leukemia (TARGET-ALL-R), and when expressed at low levels in tumor models such as uterine corpus endometrial carcinoma (TCGA-UCEC) and kidney renal clear cell carcinoma (TCGA-KIRC)(Supplementary Figure S6). ANXA1 prognosed poor outcomes when highly expressed in tumor models such as glioma (TCGA-GBMLGG), brain lower-grade glioma (TCGA-LGG), glioblastoma multiforme (TCGA-GBM), bladder urothelial carcinoma (TCGA-BLCA), uveal melanoma (TCGA-UVM), pancreatic adenocarcinoma (TCGA-PAAD), and acute lymphoblastic leukemia (TARGET-ALL-R), and when expressed at low levels in tumor models such as skin cutaneous melanoma (TCGA-SKCM) (Supplementary Figure S7). ATAD3A prognosed poor outcomes when highly expressed in tumor models such as glioma (TCGA-GBMLGG), brain lower-grade glioma (TCGA-LGG), sarcoma (TCGA-SARC), liver hepatocellular carcinoma (TCGA-LIHC), skin cutaneous melanoma (TCGA-SKCM), bladder urothelial carcinoma (TCGA-BLCA), mesothelioma (TCGA-MESO), acute myeloid leukemia (TCGA-LAML), TARGET-ALL, adrenocortical carcinoma (TCGA-ACC), and kidney chromophobe (TCGA-KICH) and when expressed in low levels in tumor models such as kidney renal papillary cell carcinoma (TCGA-KIRP) (Supplementary Figure S8).

In sum, this study established a model and identified relevant genes that can predict BLCA patient prognosis. Findings should contribute to the development of targeted BLCA treatments and help inform more precise treatment decisions.

This study is the first to assess the correlation between necroptosis-related genes and BLCA patient prognosis. Two molecular isoforms were identified by consensus clustering, immune microenvironment analysis was performed to measure prognostic ability, and three genes were identified by Lasso regression analysis. In addition, a pan-cancer analysis of the three genes was conducted to assess their prognostic relevance to other tumors. This study provides theoretical support for the development of targeted drugs and individualized therapy as well as stratified cancer treatment and evaluation.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repositories and accession numbers can be found in the supplementary material.

Acknowledgments

We thank Zheng-rong Huangand Zhong-gui Shan for technology supporting and composition instruction

Funding

This study was supported by the Natural Science Foundation of China (Grant Nos.82070291). The funder had no role in the decision to publish or preparation of the manuscript.

Availability of data and materials

Not applicable.

About this supplement

The full contents of the supplement are available online.

Author Contributions

All authors have read and approved the manuscript. Yj-L and CY developed the initial idea for this study, reviewed the completed manuscript, and supervised the process, Yl-L performed the data collection and analysis and Z-hZ prepared the manuscript, tables, and graphs. J-jX, R-jL, Z-hZ, M-wC, CY, S-xL, Y-lL, X-qW, RL, PL, and X-hZ analyzed the data. Z-gS designed the study and revised the manuscript.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations

Bray, F., Ferlay, J., Soerjomataram, I., Siegel, R.L., Torre, L.A., and Jemal, A. (2018). Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin 68, 394–424.doi.org/10.3322/caac.21492.
Chen, F., Fan, Y., Cao, P., Liu, B., Hou, J., Zhang, B., and Tan, K. (2021). Pan-Cancer Analysis of the Prognostic and Immunological Role of HSF1: A Potential Target for Survival and Immunotherapy. Oxid Med Cell Longev 2021, 5551036.doi.org/10.1155/2021/5551036.
Della Torre, L., Nebbioso, A., Stunnenberg, H.G., Martens, J.H.A., Carafa, V., and Altucci, L. (2021). The Role of Necroptosis: Biological Relevance and Its Involvement in Cancer. Cancers (Basel) 13.doi.org/10.3390/cancers13040684.
Feng, J., Lu, S.S., Xiao, T., Huang, W., Yi, H., Zhu, W., Fan, S., Feng, X.P., Li, J.Y., Yu, Z.Z., Gao, S., Nie, G.H., Tang, Y.Y., and Xiao, Z.Q. (2020). ANXA1 Binds and Stabilizes EphA2 to Promote Nasopharyngeal Carcinoma Growth and Metastasis. Cancer Res 80, 4386–4398.doi.org/10.1158/0008-5472.Can-20-0560.
Feng, Y., Feng, Q., Xu, L., Jiang, Y., Meng, F., and Shu, X. (2019). A Novel Biomarker Based on miRNA to Predict the Prognosis of Muscle-Invasive Bladder Urothelial Carcinoma. J Oncol 2019, 2654296.doi.org/10.1155/2019/2654296.
Furuya, H., Chan, O.T.M., Hokutan, K., Tsukikawa, Y., Chee, K., Kozai, L., Chan, K.S., Dai, Y., Wong, R.S., and Rosser, C.J. (2019). Prognostic Significance of Lymphocyte Infiltration and a Stromal Immunostaining of a Bladder Cancer Associated Diagnostic Panel in Urothelial Carcinoma. Diagnostics (Basel) 10.doi.org/10.3390/diagnostics10010014.
Giridhar, K.V., and Kohli, M. (2017). Management of Muscle-Invasive Urothelial Cancer and the Emerging Role of Immunotherapy in Advanced Urothelial Cancer. Mayo Clin Proc 92, 1564–1582.doi.org/10.1016/j.mayocp.2017.07.010.
Gong, Y., Fan, Z., Luo, G., Yang, C., Huang, Q., Fan, K., Cheng, H., Jin, K., Ni, Q., Yu, X., and Liu, C. (2019). The role of necroptosis in cancer biology and therapy. Mol Cancer 18, 100.doi.org/10.1186/s12943-019-1029-8.
Green, D.A., Rink, M., Xylinas, E., Matin, S.F., Stenzl, A., Roupret, M., Karakiewicz, P.I., Scherr, D.S., and Shariat, S.F. (2013). Urothelial carcinoma of the bladder and the upper tract: disparate twins. J Urol 189, 1214–1221.doi.org/10.1016/j.juro.2012.05.079.
Jin, K., Yu, Y., Zeng, H., Liu, Z., You, R., Zhang, H., Liu, C., Su, X., Yan, S., Chang, Y., Xu, L., Xu, J., Zhu, Y., and Wang, Z. (2022). CD103(+)CD8(+) tissue-resident memory T cell infiltration predicts clinical outcome and adjuvant therapeutic benefit in muscle-invasive bladder cancer. Br J Cancer 126, 1581–1588.doi.org/10.1038/s41416-022-01725-6.
Lang, L., Loveless, R., Dou, J., Lam, T., Chen, A., Wang, F., Sun, L., Juarez, J., Qin, Z.S., Saba, N.F., Shay, C., and Teng, Y. (2022). ATAD3A mediates activation of RAS-independent mitochondrial ERK1/2 signaling, favoring head and neck cancer development. J Exp Clin Cancer Res 41, 43.doi.org/10.1186/s13046-022-02274-9.
Li, F., Teng, H., Liu, M., Liu, B., Zhang, D., Xu, Z., Wang, Y., and Zhou, H. (2020). Prognostic Value of Immune-Related Genes in the Tumor Microenvironment of Bladder Cancer. Front Oncol 10, 1302.doi.org/10.3389/fonc.2020.01302.
Liu, L., Fan, J., Ai, G., Liu, J., Luo, N., Li, C., and Cheng, Z. (2019). Berberine in combination with cisplatin induces necroptosis and apoptosis in ovarian cancer cells. Biol Res 52, 37.doi.org/10.1186/s40659-019-0243-6.
Novizio, N., Belvedere, R., Pessolano, E., Morello, S., Tosco, A., Campiglia, P., Filippelli, A., and Petrella, A. (2021). ANXA1 Contained in EVs Regulates Macrophage Polarization in Tumor Microenvironment and Promotes Pancreatic Cancer Progression and Metastasis. Int J Mol Sci 22.doi.org/10.3390/ijms222011018.
Patel, V.G., Oh, W.K., and Galsky, M.D. (2020). Treatment of muscle-invasive and advanced bladder cancer in 2020. CA Cancer J Clin 70, 404–423.doi.org/10.3322/caac.21631.
Pitt, J.M., Marabelle, A., Eggermont, A., Soria, J.C., Kroemer, G., and Zitvogel, L. (2016). Targeting the tumor microenvironment: removing obstruction to anticancer immune responses and immunotherapy. Ann Oncol 27, 1482–1492.doi.org/10.1093/annonc/mdw168.
Rooney, M.S., Shukla, S.A., Wu, C.J., Getz, G., and Hacohen, N. (2015). Molecular and genetic properties of tumors associated with local immune cytolytic activity. Cell 160, 48–61.doi.org/10.1016/j.cell.2014.12.033.
Simon, N., Friedman, J., Hastie, T., and Tibshirani, R. (2011). Regularization Paths for Cox's Proportional Hazards Model via Coordinate Descent. J Stat Softw 39, 1–13.doi.org/10.18637/jss.v039.i05.
Smith, Z.L., Christodouleas, J.P., Keefe, S.M., Malkowicz, S.B., and Guzzo, T.J. (2013). Bladder preservation in the treatment of muscle-invasive bladder cancer (MIBC): a review of the literature and a practical approach to therapy. BJU Int 112, 13–25.doi.org/10.1111/j.1464-410X.2012.11762.x.
Song, Y., Liu, G., Liu, S., Chen, R., Wang, N., Liu, Z., Zhang, X., Xiao, Z., and Liu, L. (2019). Helicobacter pylori upregulates TRPC6 via Wnt/β-catenin signaling to promote gastric cancer migration and invasion. Onco Targets Ther 12, 5269–5279.doi.org/10.2147/ott.S201025.
Teng, Y., Lang, L., and Shay, C. (2019). ATAD3A on the Path to Cancer. Adv Exp Med Biol 1134, 259–269.doi.org/10.1007/978-3-030-12668-1_14.
Thakur, B., Kumar, Y., and Bhatia, A. (2019). Programmed necrosis and its role in management of breast cancer. Pathol Res Pract 215, 152652.doi.org/10.1016/j.prp.2019.152652.
Tibshirani, R. (1997). The lasso method for variable selection in the Cox model. Stat Med 16, 385–395.doi.org/10.1002/(sici)1097-0258(19970228)16:4<385::aid-sim380>3.0.co;2-3.
Walter, W., Sánchez-Cabo, F., and Ricote, M. (2015). GOplot: an R package for visually combining expression data with functional analysis. Bioinformatics 31, 2912–2914.doi.org/10.1093/bioinformatics/btv300.
Wang, X., Bai, X., Yan, Z., Guo, X., and Zhang, Y. (2021). The lncRNA TUG1 promotes cell growth and migration in colorectal cancer via the TUG1-miR-145-5p-TRPC6 pathway. Biochem Cell Biol 99, 249–260.doi.org/10.1139/bcb-2020-0017.
Wilkerson, M.D., and Hayes, D.N. (2010). ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics 26, 1572–1573.doi.org/10.1093/bioinformatics/btq170.
Witjes, J.A., Compérat, E., Cowan, N.C., De Santis, M., Gakis, G., Lebret, T., Ribal, M.J., Van Der Heijden, A.G., and Sherif, A. (2014). EAU guidelines on muscle-invasive and metastatic bladder cancer: summary of the 2013 guidelines. Eur Urol 65, 778–792.doi.org/10.1016/j.eururo.2013.11.046.
Yang, G., Zhang, L., Liu, M., Liu, Q., Duan, X., and Bo, J. (2019). CD163 + macrophages predict a poor prognosis in patients with primary T1 high-grade urothelial carcinoma of the bladder. World J Urol 37, 2721–2726.doi.org/10.1007/s00345-018-02618-1.
Yang, L.L., Liu, B.C., Lu, X.Y., Yan, Y., Zhai, Y.J., Bao, Q., Doetsch, P.W., Deng, X., Thai, T.L., Alli, A.A., Eaton, D.C., Shen, B.Z., and Ma, H.P. (2017). Inhibition of TRPC6 reduces non-small cell lung cancer cell proliferation and invasion. Oncotarget 8, 5123–5134.doi.org/10.18632/oncotarget.14034.
Yoshihara, K., Shahmoradgoli, M., Martínez, E., Vegesna, R., Kim, H., Torres-Garcia, W., Treviño, V., Shen, H., Laird, P.W., Levine, D.A., Carter, S.L., Getz, G., Stemke-Hale, K., Mills, G.B., and Verhaak, R.G. (2013). Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun 4, 2612.doi.org/10.1038/ncomms3612.
Yu, G., Wang, L.G., Han, Y., and He, Q.Y. (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. Omics 16, 284–287.doi.org/10.1089/omi.2011.0118.
Zhao, Y., Li, Y., Wang, L., Yang, H., Wang, Q., Qi, H., Li, S., Zhou, P., Liang, P., Wang, Q., and Li, X. (2013). microRNA response elements-regulated TRAIL expression shows specific survival-suppressing activity on bladder cancer. J Exp Clin Cancer Res 32, 10.doi.org/10.1186/1756-9966-32-10.

No competing interests reported.

Download PDF

Version 1

posted

You are reading this latest preprint version

Using Necroptosis-Associated Genes To Predict The Immune Microenvironment And Prognosis Of Bladder Urothelial Carcinoma

Status:

Version 1

Abstract

PURPOSE

METHODS

RESULTS

CONCLUSION

Figures

1 Introduction

2. Materials And Methods

2.1 Patients and Datasets

2.2 Screening for prognosis-related genes

2.3 Cluster analysis

2.4 TIME evaluation and immune infiltration analysis

2.5 Functional enrichment analysis

2.6 Risk model building and validation

2.7 Pan-cancer analysis

2.8 Data analysis

3. Results

3.1 Clustering by necroptosis-related genes identified two molecular subtypes

3.2 Immunological analysis

3.3 DEG expression and functional analyses

3.4 Training queue risk modeling

3.5 Screening for independent predictors of BLCA patient prognosis

3.6 BLCA risk model relationship to the TIME and disease prognosis

3.7 Construction and calibration of a BLCA integrated model

4. Discussion

5. Conclusion

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1