ScRNA-seq and bulk RNA-seq reveal the characteristics of NK cell and establish a risk signature in bladder cancer

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

Download PDF

Research Article

ScRNA-seq and bulk RNA-seq reveal the characteristics of NK cell and establish a risk signature in bladder cancer

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

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Bladder cancer poses a significant risk of recurrence and progression to intramuscular disease, also known as non-muscle invasive bladder cancer (NMIBC). Immunotherapy has emerged as a pioneering strategy for managing non-muscle invasive bladder cancer (NMIBC). However, there exist variations in the resistance patterns of immune checkpoint inhibitors (ICIs) and the corresponding patient response, thereby necessitating further investigation into novel therapeutic interventions.Several research studies have indicated that the focused metabolism of sulphur dioxide could potentially serve as a substitute for existing approaches to treating bladder cancer. The research was centred on the examination of genetic variations of genes responsible for the metabolism of sulphur dioxide (DRGs) that are linked to bladder cancer. The research revealed that out of the 412 samples analysed, 136 exhibited mutations in DRGs, indicating a mutation rate of 33.01%.The research additionally discovered that duplications or deletions of genomic segments, known as copy number variants (CNV), are prevalent in bladder carcinoma, particularly in relation to dorsal root ganglia (DRGs).The investigation additionally examined the impact of genetic variations on messenger RNA (mRNA) expression. The results indicated a positive correlation between the level of expression of specific genes in bladder cancer tissue and the variation in replicating numbers. The investigation also recognised two types of sulphur dioxide metabolism that are dispersed, and are linked with distinct gene expression patterns and characteristics of immune cell infiltration. The investigation enhanced the examination of genes with differential expression and revealed that the bladder cancer tumour microenvironment was linked to the metabolism of sulphur dioxide, which resulted in varying prognostic outcomes. The investigation additionally formulated a prognostic instrument utilising clinical parameters and risk assessments. The potential of the predictive signature to correlate with the metabolism of sulphur dioxide has implications for the characterization of a patient's immune landscape and the prediction of their prognosis. However, additional research is required to address various concerns, such as the validation of models' accuracy and efficacy, as well as the undertaking of more prospective and fundamental investigations.

Sc and buck RNA-seq analysis

risk score signature

NK cell

BLCA

Bladder cancer, also known as BLCA, ranks as the sixth most prevalent form of cancer in males. This type of cancer is associated with significant risks of both recurrence and mortality[1]. According to a report, it was estimated that in the United States, the incidence of new cases of BLCA would be approximately 81,400, and the mortality rate would be approximately 17,980 in the year 2020[2]. At present, cisplatin-based chemotherapy is the primary treatment modality for advanced BLCA[3, 4]. However, the response duration to chemotherapy is restricted and the prognosis for patients is unfavourable, with a 5-year relative overall survival rate of only 15% for metastatic disease[5]. According to Cathomas et al.[6], in cases where the disease remains localised, recommended treatment options for these patient groups include radical cystectomy with urinary diversion, radical radiotherapy, and neoadjuvant and/or adjuvant chemotherapy or immunotherapy. According to Powles et al.[7], individuals with advanced metastatic UCB may experience improved overall survival (OS) through cisplatin-containing chemotherapy. Additionally, initial results indicate that immunotherapy utilising PD-1/PD-L1 checkpoint inhibitors may be a promising treatment option. Nevertheless, these interventions entail notable surgical hazards and toxicities associated with the treatment. Despite the advancements in therapeutics, patients diagnosed with advanced or invasive bladder cancer (BLCA) continue to experience significant morbidity and unfavourable clinical outcomes. At present, there exists a dearth of tools that can aid in the identification of patients who are more likely to derive therapeutic benefits from chemotherapy or immunotherapy within these cohorts. Hence, there is a pressing need for novel biomarkers that can function as therapeutic targets or prognostic indicators to aid in distinguishing between individuals who respond to treatment and those who do not.

According to a recent research, cells with high levels of SLC7A11 may prevent NK cell when glucose is scarce through cystine absorption, which is mediated by SLC7A11[8, 9]. Nonetheless, this phenomenon may potentially give rise to a novel form of cellular demise known as NK cell. Recent research indicates that NK cell can induce apoptosis in tumour cells by modulating the conformation of cytoskeletal proteins and is associated with modifications in the cellular redox status. Consequently, NK cell has the potential to introduce a novel avenue in the management of neoplastic growths. However, additional research and analysis are imperative to gain a comprehensive understanding of its distinctive mechanisms and clinical applications. There has been a notable association established between cancer and disulfide metabolism in recent times. Disulfide metabolism pertains to the intracellular redox mechanisms primarily regulated by the formation and cleavage of disulfide bonds. According to recent research, cancer cells frequently experience oxidative stress, leading to disulfide metabolism complications that affect their survival and proliferation [10, 11]. The metabolic process of disulfide in cancer cells has been found to be associated with various biological phenomena such as immunological evasion, metastasis, and resistance to treatment[12, 13]. Disulfidoptosis, a newly identified form of programmed cell death, may be associated with the immune response to tumours. The signal of cell death caused by NK cell has the potential to activate tumour immune cells, which can subsequently stimulate the immunological response of tumor-specific NK cells. This process can enhance both humoral and cellular immunity, ultimately leading to improved effectiveness of cancer treatments. It is imperative to identify novel biomarkers linked to NK cell and establish additional associations between disulfide metabolism-related targets and pathways with cancer susceptibility.

The present investigation examined the genetic variations in genes associated with NK cell (DRGs) in a cohort of individuals diagnosed with bladder cancer. Among the 412 specimens that were analysed, 136 exhibited mutations in DRGs, resulting in a mutation frequency of 33.01%. The gene MYH10 exhibited the highest frequency of mutations, while MYH9, FLNB, and ACTB followed suit in descending order. Copy number variations (CNVs) are commonly observed in bladder cancer, with a considerable proportion of differentially regulated genes (DRGs) exhibiting either copy number amplifications or deletions. The research investigation additionally examined the influence of genetic variations on messenger RNA expression in the context of bladder cancer. The findings of the study revealed that copy number variation could potentially serve as a mediator for variations in dorsal root ganglion expression. The study identified two discrete phenotypes of disulfidoptosis, which were correlated with specific gene expression patterns and immune cell infiltration characteristics. The correlation between patient prognosis and the presence of NK cell clusters A and B was observed, where cluster A exhibited a positive association with improved outcomes. Furthermore, a secondary clustering approach was employed based on the expression patterns of genes that were differentially expressed. This analysis identified three distinct gene clusters that were found to be associated with NK cell phenotypes. The predictive ability of patient prognosis in various cohorts was determined by computing a NK cell score using gene expression levels. The assessment comprised of a set of ten genes that were classified into either high-risk or low-risk categories. The score for NK cell exhibited noteworthy variations among various modes of clustering and proved to be effective in forecasting survival results in the datasets for training, complete, and testing. A nomogram was utilised to integrate the score, thereby enhancing the prognostic accuracy for individuals diagnosed with bladder cancer. The investigation also examined the correlations between the NK cell score and tumour mutation burden, along with genomic instability, underscoring the importance of genetic mutations in the progression of tumours.

Bladder Cancer Datasets and Preprocessing

The current study involved an assemblage of 1020 specimens of bladder cancer (BLCA) obtained from five high-throughput platforms. The investigation employed a cumulative sum of 1,020 specimens procured from diverse databases, comprising 412 specimens from The Cancer Genome Atlas (TCGA) database, 189 specimens from GSE13507, 93 specimens from GSE31684, 308 specimens from GSE32894, and 18 specimens from GSE37815. Furthermore, a heterogeneous array of data was collected, comprising of somatic mutation data, copy number variation (CNV) data, clinical particulars, and survival data. All study participants had complete follow-up information and available RNA-sequencing (RNA-seq) data. Participants who did not have complete survival data were excluded from subsequent analysis. The study employed normalised matrix files sourced from GEO and gene expression data measured in Fragments Per Kilobase Million (FPKM values) obtained from TCGA for analysis purposes. The transformation of FPKM values into TPM values was carried out to facilitate integrated analysis. To address batch effects and non-biological deviations, we employed the "ComBat" algorithm[14] which is implemented in the R package "SVA". The single-cell sequencing data of seven untreated muscle-infiltrating bladder cancers from GSE135337[15]was also obtained. The data underwent analysis in accordance with the prescribed post-quality control procedures, and the cells were annotated based on the information presented in the original study. The functions "FindAllMarkers" and "FindMarkers" are utilised to execute Wilcoxon tests between cell cluster pairs with the aim of identifying genes that are expressed in a cluster-specific manner. The function 'featureplot' is utilised to visually represent the expression of a particular gene.

Unsupervised Clustering Based on NK cell-Related Genes.

In order to discern unique NK cell patterns linked to NK cell-related genes (DRGs), a comprehensive collection of 26 DRGs was collated from previously published scholarly works. The samples underwent hierarchical agglomerative clustering using the R package "ConsensusClusterPlus"[16]in a subsequent step. The aforementioned methodology utilised stability indicators within unsupervised analysis to ascertain the most suitable quantity of clusters and their respective constituents. In order to guarantee the robustness of the clustering outcomes, the procedure was repeated for 1,000 iterations, thus yielding dependable stability metrics.

Gene Set Variation Analysis

The R package "GSVA" [17]was utilised to perform enrichment analysis aimed at exploring the differences in biological processes across the various subtypes of NK cell. The utilisation of the gene set variation analysis (GSVA), which is a non-parametric and unsupervised approach, was implemented to evaluate the fluctuations in enrichment pertaining to pathways and biological process activities among diverse expression datasets. The GSVA analysis employed Gene Ontology (GO) and Kyoto Encyclopaedia of Genes and Genomes (KEGG) gene sets sourced from the MSigDB database (http://software.broadinstitute.org/gsea/msigdb/). The NK cell-related pathways were visualised using heatmaps, and statistical significance was established by applying an adjusted p-value threshold of less than 0.05.

Estimating of Immune Cell Infiltration

In order to evaluate and measure the extent of immune cell infiltration (ICI) in each specimen, the methodology of single-sample gene-set enrichment analysis (ssGSEA) was utilised, with the aid of the R package "GSVA". Initially, immune cell marker gene expression data was acquired from Charoentong's investigation. The ssGSEA technique was employed to compute the enrichment score, which functioned as a representative indicator of the proportional prevalence of immune cell infiltration in every specimen. Through the comparison of these scores, we conducted an analysis of the ICI associated with the different subtypes of NK cell.

Identification of Differentially Expressed Genes Between NK cell Subtypes in Bladder Cancer

The R package "Limma"[18]was employed to identify differentially expressed genes (DEGs) among various NK cell subtypes in bladder cancer. Differentially expressed genes (DEGs) were identified using a significance threshold of adjusted p-values less than 0.01. In addition, the authors performed Gene Ontology (GO) and Kyoto Encyclopaedia of Genes and Genomes (KEGG) functional annotation analyses using the R package "clusterProfiler"[19] to enhance their comprehension of the functional implications linked to the differentially expressed genes (DEGs).

Construction of NK cell Score

In order to evaluate the predictive significance of DEGs, a univariate COX regression analysis was employed to examine the DEGs that are linked to prognosis. The individuals under study were categorised into gene clusters A and B for further examination. The patients were randomly allocated into two groups, namely a training set (n = 452) and a testing set (n = 452), in a 1:1 ratio. The risk score linked with DEGs was developed using the training set. The LASSO method was employed to address the issue of overfitting and to reduce the number of differentially expressed genes (DEGs) under consideration, ultimately resulting in the identification of a set of 10 genes suitable for the development of a prognostic model. The risk score, also known as the NK cell score, was computed utilising the subsequent mathematical expression: Risk score = Σ (Expi * Coefi), where Coefi represents the risk coefficient and Expi denotes the gene expression. Following the computation of risk, the patients were stratified into two groups, namely low-risk and high-risk, based on the median risk score. Subsequently, survival analysis and receiver operating characteristic (ROC) curves were conducted on the training set. Validation was conducted on both the testing and entire sets. The nomogram was built using the “RMS” package of R software to predict individual survival probability, and calibration curves for the prediction of 1 -, 3-, and 5-year survival rate of BLCA patients were plotted.

Prediction of Multiple Therapeutic Sensitivities

The study aimed to analyze the differences in therapeutic sensitivities between the high risk group and low risk group concerning targeted inhibitor (TI) therapy and immunotherapy. The concentration causing a 50% reduction in growth (IC50) of TIs, including Notch, Hedgehog (HH), and Wnt inhibitors, was calculated using the R package "pRRophetic"[20]. The Wilcoxon rank-sum test was employed to compare the differences in IC50 values between high risk group and low risk group. For evaluating immunogenicity, multiple genes associated with effector cells, MHC molecules, immune regulatory factors, and immunosuppressive cells were considered. Immunogenicity was estimated and quantified using machine learning methods to minimize bias. The IPS values were compared between high risk group and low risk group using different immunotherapy methods to predict immunotherapy sensitivity. Furthermore, publicly available datasets treated with immunotherapy were comprehensively searched. The metastatic urothelial tumors cohort (IMvigor210:http://research-pub.gene.com/IMvigor210CoreBiologies)[21]. The IMvigor210 data were preprocessed using the R package "IMvigor210CoreBiologies." The RNA-seq data were filtered and normalized using the R package "edgeR"[22], and the data were transformed using the "voom" function in the R package "limma". Prognostic information and immunotherapy outcomes were also collected. The NK cell Scores of each sample in the IMvigor210 cohort was also quantified.

Analysis of ten hub NK cell score related genes by scRNA-seq

An analysis of single-cell RNA sequencing data was performed using R packages, namely "Seurat" and "SingleR". In order to ensure the preservation of scRNA-seq data of superior quality, the raw matrix for each cell underwent three filtering measures. These measures entailed the inclusion of solely those genes that were expressed in a minimum of 5 single cells, the elimination of cells that expressed fewer than 100 genes, and the removal of cells that exhibited over 5% of mitochondrial genes. The scRNA-seq data was initially normalised using the "Seurat" R package, wherein the "NormalizeData" function was employed with the normalisation method being set as "LogNormalize." The scRNA-seq data was normalised and converted into a Seurat object. The "FindVariableFeatures" function was utilised to identify the top 1,500 genes with high variability. Subsequently, we utilised the "RunPCA" function from the "Seurat" R package to execute principal component analysis (PCA) on the scRNA-seq data, with the aim of reducing its dimensionality. This was based on the top 1,500 genes. The JackStraw method was employed to identify principal components (PCs) that were statistically significant. Based on the proportion of variance explained, we selected the initial 15 PCs for the purpose of conducting cell clustering analysis. Cell clustering analysis was conducted using the "Seurat" package's "FindNeighbors" and "FindClusters" functions. The k-nearest neighbour graph was generated by utilising the Euclidean distance metric in principal component analysis (PCA). The "FindNeighbors" function was employed to identify the nearest neighbours of each cell. Subsequently, the "RunTSNE" function was utilised to execute t-distributed stochastic neighbour embedding (t-SNE). The utilisation of t-SNE-1 and t-SNE-2 facilitated the demonstration of cell clustering. The Wilcoxon-Mann-Whitney test was utilised to compute the differentially expressed genes (DEGs) of each cluster through the "FindAllMarkers" function present in the "Seurat" package. The marker genes for each cluster were identified using a cutoff threshold of adjusted p-value < 0.01 and |log2 (fold change)| >1. The process of cluster annotation involved the utilisation of reference data from the Human Primary Cell Atlas for reference-based annotation.

Statistical Analyses

The statistical analyses conducted in this study were executed utilising the R statistical language (version 4.0.5). The Wilcoxon test was employed to conduct a comparison between the two groups. The Kruskal-Wallis test was employed to conduct comparisons among multiple groups. The Kaplan-Meier plotting tool was utilised to generate prognostic survival curves, while the log-rank test was employed to evaluate the statistical significance of any observed differences. The Spearman correlation test was employed to conduct correlation analysis and determine the correlation coefficient. The R package "maftool"[23] was utilised to demonstrate the gene mutation status across various groups. The R package "rcircos"[24] was utilised to cartographically represent the copy number variation (CNV) of DRGs across 23 pairs of chromosomes. A significance level of p < 0.05 was utilised in all analyses.

Genetic Variations in NK cell-Related Genes among bladder cancer

The current study examined a group of 26 genes that are linked to NK cell, which are commonly known as NK cell-Related Genes (DRGs), based on recent academic literature. The examination of 412 specimens indicated that 136 of them exhibited mutations in DRGs, leading to a mutation frequency of 33.01%. The findings suggest that MYH10 exhibited the most frequent incidence of mutations, trailed by MYH9, FLNB, and ACTB, as depicted in Fig. 1A. The results of our analysis suggest that Copy Number Variation (CNV) aberrations are frequently detected in bladder cancer within the context of Disease-Related Genes (DRGs). The findings depicted in Fig. 1B indicate that a significant fraction of DRGs manifested copy number deletions, while nearly half of the DRGs demonstrated a higher incidence of CNV amplification. The results mentioned above emphasise the significant occurrence of copy number variations (CNV) and the unique characteristics of deletions or amplifications in various disease-related genes (DRGs). It is worth mentioning that certain genes exhibited a higher frequency of copy number amplification, namely SLC3A2, ACTN4, OXSM, and TLN1. By way of comparison, it was noted that specific genes, specifically FLNB, NDUFS1, PDLIM1, and NDUFS1, displayed a higher incidence of deletions. Furthermore, a study was carried out to ascertain the incidence of copy number variations (CNV) within the dorsal root ganglia (DRGs) situated on the chromosome, as illustrated in Fig. 1C. A study was conducted to assess the impact of genetic variations on mRNA expression in BLCA. This involved an analysis of mRNA expression levels in both healthy and malignant tissues. The results suggest that copy number variation (CNV) may act as a mediator for variations in the expression of dorsal root ganglion (DRG). The findings suggest that DRGs with a higher frequency of amplification demonstrate a significant increase in expression within cancerous tissues, encompassing SLC3A2, OXSM, and GYS1, among others. By contrast, it was observed that certain genes, namely NDUFS1 and NCKAP1, exhibited a decrease in their expression levels in malignant tissues relative to healthy tissues, as illustrated in Fig. 1D. The observed heterogeneity in the expression of DRGs in both normal and malignant tissues underscores the crucial importance of distinct DRGs expression in the initiation and advancement of BLCA.

Two Distinct NK cell Phenotypes Associated with NK cell-Related Gene Expression Patterns and Immune Cell Infiltration Characteristics and Biological Behaviours

Univariate Cox proportional hazards regression and correlation analyses were performed to gain further understanding of the interrelationships, associations, and prognostic significance of DRGs. The results of the study suggest that the prognosis of patients with bladder cancer (BLCA) was significantly influenced by 18 functional-related genes (FRGs), all of which had p-values below 0.05. It is worth mentioning that a positive correlation was identified between DRGs that displayed a similar prognostic impact. The findings suggest that there exists a noteworthy and favourable association between the manifestation of FLNA and the manifestation of MYH9, MYH10, ACTN4, TLN1, MYL6, IQGAP1, ACTB, DSTN, and CAPZB. Conversely, a significant negative correlation was observed among functional response groups (FRGs) that displayed favourable and unfavourable prognostic factors. The advantageous factor OXSM exhibited an inverse correlation with the risk prognostic factors ACTB, TLN1, and FLNA. The results suggest that TLN1, a factor associated with increased risk, demonstrated a negative correlation with the advantageous prognostic factors CD2AP, NDUFA11, and OXSM, as depicted in Fig. 2A. The findings indicate the presence of complex interconnections between DRGs, which carry important implications for patient prognosis and the detection of TME cell infiltration. Cluster analysis was performed on the distinct DRG expression patterns of patients using the R package "ConsensusClusterPlus". The current investigation has effectively recognised two separate phenotypes associated with NK cell, denoted as NK cell clusters A and B (as depicted in Fig. 2B and 2C). The findings of the prognostic evaluation revealed that patients classified as NK cell cluster A exhibited a significantly better prognosis as compared to those assigned to NK cell cluster B (p < 0.001) (Fig. 2D). Furthermore, a heatmap was employed to visually represent the expression profile of dorsal root ganglia (DRGs), highlighting their differential expression in the two clusters (refer to Fig. 2E). A Kyoto Encyclopaedia of Genes and Genomes (KEGG) Gene Set Variation Analysis (GSVA) was conducted to examine the differences in biological behaviours among the NK cell clusters. The results suggest that there was a significant increase in stromal-related signals within the NK cell cluster B. The signals mentioned earlier were composed of the MAPK signalling pathway and ECM receptor interaction, which are illustrated in Fig. 2F. The findings of the principal component analysis (PCA) demonstrate that genes associated with NK cell have the ability to differentiate between patients with bladder cancer (BLCA), as depicted in Fig. 2G. The analysis of immune cell infiltration (ICI) revealed that the NK cell cluster B exhibited a significant presence of various immune cells, including B cells, CD4 + NK cells, CD8 + NK cells, Eosinophils, MDSCs, Macrophages, MasNK cells, Neutrophil, and regulatory NK cells (as depicted in Fig. 2H). The results mentioned above are noteworthy as they indicate that disulfidoptosis may have the capability to elicit immunosuppressive and inflammatory reactions within the tumour microenvironment (TME).

Secondary Clustering Using Differentially Expressed Genes and Identification of Prognostic-Related Subtypes

To perform a more thorough examination of the potential biological behaviours associated with each NK cell phenotype, a set of 1213 differentially expressed genes (DEGs) was detected in relation to the NK cell phenotypes using the R package "Limma" (Fig. 3A). The differentially expressed genes (DEGs) were subjected to enrichment analyses utilising GO and KEGG. The findings suggest a noteworthy enhancement in diverse immune-associated biological pathways, thus corroborating the pivotal implication of NK cell in the tumour microenvironment (TME). The aforementioned results are depicted in Figs. 3B and 3C. Univariate Cox analyses were conducted to assess the prognostic relevance of 1213 genes that exhibited differential expression. In the end, a comprehensive set of 923 genes exhibiting differential expression (DEGs) were pinpointed and chosen for further scrutiny, in relation to their prognostic implications for patients afflicted with bladder cancer. To bolster the legitimacy of this phenomenon, an unsupervised cluster analysis was conducted using the 923 differentially expressed genes (DEGs). The outcome of this study led to the categorization of patients into three distinct clusters (gene clusters A and B), which exhibited a correlation with NK cell phenotypes (as depicted in Fig. 3D and 3E). The stratification was implemented with the aim of capturing additional phenotypic variations that result from disulfidoptosis patterns. The findings of this study demonstrate that the utilisation of the clustering algorithm demonstrated that patients categorised under gene cluster A displayed a significantly favourable prognosis (p < 0.001) (Fig. 3F). The results indicate that a significant proportion of patients in gene cluster A were linked to NK cell cluster A, while those in gene clusters B, which exhibited a relatively unfavourable prognosis, were linked to NK cell cluster B (as shown in Fig. 3G). The investigation unveiled significant discrepancies in the presentation of DRGs among discrete gene clusters, as illustrated in Fig. 3H. The results align with the expected effects associated with disulfidoptosis patterns.

The NK cell Score and its Predictive Power for Prognosis in Multiple Cohorts

The score for NK cell was computed through the utilisation of differential expression of genes (DEGs) that are associated with NK cell. The cohort was subjected to a random partitioning procedure, which led to the creation of two distinct sets, namely the training set (n = 452) and the testing set (n = 452). Subsequently, a LASSO regression analysis was conducted on a cohort of 923 genes that were previously identified as exhibiting differential expression (DEGs) to ascertain the optimal prognostic signatures. The utilisation of LASSO regression yielded the discovery of ten genes that exhibit a correlation with NK cell, as depicted in Figs. 4A and 4B. Subsequently, a multivariate Cox regression analysis was conducted on the 10 genes that have been associated with NK cell. The findings of the analysis yielded a definitive classification of ten genes, namely PPP1R3C, SLC12A8, SCARA3, NR2F1, SLC7A5, DUSP2, SLCO1B3, CD3D, OLFM4, and UGT1A1. The dataset under consideration comprises of ten genes, out of which six have been categorised as high-risk genes, namely SLC12A8, SCARA3, NR2F1, SLC7A5, SLCO1B3, and OLFM4. The remaining four genes have been classified as low-risk genes, namely PPP1R3C, DUSP2, CD3D, and UGT1A1. The approach employed for computing the risk score was as follows: The process of calculating the risk score involves the evaluation of the levels of expression of several genes, namely PPP1R3C, SLC12A8, SCARA3, NR2F1, SLC7A5, DUSP2, SLCO1B3, CD3D, OLFM4, and UGT1A1. The coefficients of the gene are presented as follows: -0.101162561261356, 0.323897733590379, 0.187765237313683, 0.335582585946667, 0.248807225420312, -0.157571966578649, 0.149888375776253, -0.587488152209679, 0.090212513648201, and − 0.256232444902395. The distribution of two subtypes, two gene clusters, and two NK cell score groups is depicted in Fig. 4C. Furthermore, a research study was conducted to examine the correlation between the NK cell score and the two clustering modes, which demonstrated notable variations in the NK cell score across the different clustering modes. The statistical analysis reveals a significant difference in the NK cell score between NK cell cluster A and NK cell cluster B (p < 0.001), as illustrated in Fig. 4D. The statistical analysis results reveal that the NK cell score median of gene cluster B was significantly greater than that of gene clusters A (p < 0.001), as illustrated in Fig. 4E. The study revealed notable inconsistencies in the portrayal of Diagnosis-Related Groups (DRGs) across various NK cell scores, as depicted in Fig. 4F. The results obtained from the survival analysis indicate a noteworthy distinction between the cohorts with elevated and diminished NK cell scores within the training dataset. The statistical significance of the findings is corroborated by the data presented in Fig. 5A, where the p-value is reported to be less than 0.001. The diagrammatical depiction presented in Fig. 5B elucidates the Area Under the Curve (AUC) metrics, which serve as an indicator of the efficacy of different risk scores in prognosticating rates for intervals of 1, 3, and 5 years in the training cohort. The AUC values for the three instances, as determined by the area under the ROC curve, are 0.808, 0.844, and 0.845, respectively. The results of the survival analysis indicate that the low-risk cohort demonstrated a significantly better prognosis as compared to the high-risk cohort. This suggests that the disulfidoptosis score is a highly effective prognostic indicator in instances of bladder cancer (BLCA). To assess the prognostic predictive capacity of the risk model, we conducted a replication of the procedures in both the complete set (as presented in Figs. 5C-D) and the testing set (as demonstrated in Figs. 5E-F). The results obtained were in accordance with those obtained in the training dataset, as demonstrated by the risk score and survival plots. ROC curves were produced for both the complete cohort and the testing set in order to assess the predictive capacity of various risk scores in relation to the sensitivity and specificity of 1-, 3-, and 5-year BCR rates. The study reports the AUC values for the occurrence of at 1-, 3-, and 5-year intervals in the entire cohort and the testing set. The AUC values for the entire cohort were 0.764, 0.764, and 0.763 for the respective intervals, while the corresponding values for the testing set were 0.724, 0.692, and 0.690. The study findings indicate notable variations in survival outcomes among patients categorised into high and low risk groups, despite the absence of significant distinctions in age and gender, as illustrated in Fig. 5G. A nomogram was developed with the objective of providing clinicians with an improved tool for the quantitative prognostication of patients with BLCA. This was achieved by combining various factors, including grade, M, gender, risk scores, age, N, stage, and T. The nomogram, a visual aid, demonstrated the significant relevance of the risk score with respect to various clinical parameters, as depicted in Fig. 5K. Furthermore, calibration curves were constructed to evaluate the correspondence between the anticipated and observed survival results of individuals diagnosed with BLCA, demonstrating a positive concurrence between the nomogram and the factual survival percentages (Fig. 5L). The findings of the study indicate that the Nomogram exhibited superior net benefit for 1-year and 3-year overall survival (OS), as well as for 5-year OS, as demonstrated by the results of the decision curve analysis (DCA) (refer to Fig. 5H-6J). This study provides evidence that the nomogram, which integrates the risk scores derived from our analysis, displays a notable level of precision in predicting the overall survival of individuals who have been diagnosed with bladder cancer.

The NK cell Score and its Associations with Tumor Mutation Burden and Genomic Instability

The significance of genetic mutations in the development of tumours is widely recognised in the field of oncology. The current study entailed an examination of somatic mutation data utilising disulfidoptosis score obtained from the TCGA database. The results of the analysis revealed that TP53, TTN, and KMT2D were the three genes that were most commonly mutated in both the high risk and low risk groups (as depicted in Fig. 6A and 6B). Furthermore, the subjects were categorised into high and low cohorts based on the median TMB value, followed by the execution of a Kaplan-Meier survival analysis. The study's results indicate that patients categorised under the high TMB group demonstrated a more favourable prognosis, implying that TMB has the potential to function as a prognostic marker for BLCA patients, as illustrated in Fig. 6C. To perform a more thorough examination of the effect of risk score and TMB on patient survival, the patient cohort was divided into four distinct subgroups and underwent a rigorous survival evaluation. The study's results indicate that the subgroup characterised by low TMB and high risk score had the least favourable prognosis, thus confirming the model's effectiveness and identifying the most favourable prognostic subgroup for clinical use (as shown in Fig. 6D). In summary, our analysis of somatic mutations and tumour mutational burden (TMB) confirms the prognostic significance of the disulfidoptosis score and provides valuable insights into the underlying mechanisms implicated in the progression and development of BLCA.

Immune Cell Infiltration Characteristics and Biological Behaviours between the high risk and low risk group

The cibersort algorithm was utilised to estimate the composition of 22 infiltrating immune cells in both high and low risk patient groups (as depicted in Figs. 7A and 7B). Figure 7C depicts the interrelationships among the 22 immune cells. The study conducted an analysis of the differential proportion of infiltrating immune cells between groups categorised as high risk and low risk. The results depicted in Fig. 7A demonstrate that, in comparison to the low risk group, the high risk group exhibited significant downregulation of CD8 NK cells (P < 0.001) and Tregs (P < 0.01), while displaying significant upregulation of M0 macrophages (P < 0.01) and M2 macrophages (P < 0.01). The study employed Pearson correlation analysis to examine the association between ten genes (PPP1R3C, SLC12A8, SCARA3, NR2F1, SLC7A5, DUSP2, SLCO1B3, CD3D, OLFM4, and UGT1A1) linked to disulfidoptosis and 22 immune cell types, as depicted in Fig. 7D. The Violin Plot analysis revealed that there was a statistically significant difference between the high-risk and low-risk groups in terms of stromal, immune, ESTIMATE, and tumour purity scores. Specifically, the stromal, immune, and ESTIMATE scores showed significant differences between the two groups (refer to Fig. 7G). The ssGSEA method was utilised to examine the associations between the risk scores and immune-related cellular components and pathways. The findings indicate notable dissimilarities between the low-risk and high-risk cohorts (as illustrated in Figs. 7E), specifically in relation to B cells, pDCs, Th2 cells, and TIL. The analysis of high-risk and low-risk groups within the immune subtype groups revealed a statistically significant difference in the risk score (P < 0.001) (Fig. 7F). Furthermore, examination of tumour stemness revealed that RNAss exhibited a decrease in expression levels within the high-risk cohort, and a negative correlation was observed between the risk scores and RNAss (as depicted in Fig. 7H).

The NK cell Score Had Great Potential for Predicting Immunotherapy Efficacy

The present study aimed to evaluate the predictive capacity of the NK cell score in relation to the effectiveness of immunotherapy for patients with bladder cancer. The IMvigor 210 cohort of anti-PD-L1 immune therapy was utilised to forecast the reaction to immune therapy. There were notable variations in the risk scores across the response groups, wherein the LR group exhibited a greater percentage of CR/PR, as depicted in Figs. 8A and B. The area under the receiver operating characteristic curve (AUC) values obtained from the IMvigor 210 study were reported to be 0.568, as depicted in Fig. 8C. Nevertheless, the ICB therapy demonstrates efficacy solely in a specific subgroup of patients. In order to conduct a more comprehensive investigation into the function of risk scores in the context of immune therapy, we proceeded to utilise the TIDE score as a means of assessing the prospective immune dysfunction present in both tumours and regional lymph nodes. The findings indicate that individuals categorised as low risk exhibited a greater likelihood of responding to immune therapy, as demonstrated by Figs. 8D-8J. The present study aimed to investigate the correlation between the risk scores and the IC50 values of two frequently prescribed chemotherapy agents for bladder cancer. The results indicate that Vinblastine and Cisplatin exhibit greater sensitivity in the low risk group, as illustrated in Figs. 8K-8L. In general, the results indicate that the disulfidoptosis score holds considerable promise in forecasting the prognosis and effectiveness of immunotherapy among individuals with bladder cancer (BLCA). The association between the NK cell score and differential drug sensitivity, expression of immune checkpoint genes, and clinical response to immunotherapy has been established.

Analysis of scRNA-seq Data of ten hub NK cell score related genes in BLCA

We used the bladder cancer single-cell dataset GSE135337 to analyze the expression of 10 hub NK cell score related genes in the tumor microenvironment (TME). After performing gene filtering, normalization, principle component analysis, we first set the resolution at 0.7, resulting in the 37 cell clusters, and visualized the original samples utilizing t-SNE plots (Figs. 9A-9C). DEGs for cell types were identified (Fig. 9D). Based on the expression patterns of markers from CellMarker, we manually annotated these clusters as the following 5 cell types (Fig. 9E). Then, to explore the expression patterns of ten hub NK cell score related genes in BLCA at the single-cell level, we drew t-SNE plots of these hub NK cell score related genes. In general, only seven hub NK cell score related genes were expressed in all 5 cell types, but not in all cells (Figs. 9F-9H). Subsequently, we explored and predicted the interactions between 5 cell types in the TME. First, we calculated the interactions between different types of macrophages and other cells, and assessed the interaction strength (Fig. 9I). In addition to their self-interaction, the interactions between all kinds of cell types with monocyte cells were the most obvious and frequent (Fig. 9J). Bubble plot of cell-to-cell communication between different subtypes of 5 cell types was showed in Fig. 9K.

BLCA is a prevalent urological malignancy that is observed globally, with an estimated incidence of around 550,000 new cases annually [25]. Urothelial carcinomas of the bladder (UCB) account for more than 90% of bladder cancers. The evolution of UCB's management has undergone significant changes in recent decades. Non-muscle invasive bladder cancer (NMIBC) is a prevalent condition among patients, accounting for around 70% of cases. The standard treatment for this type of cancer involves transurethral resection of bladder tumour (TURBT), followed by intravesical chemotherapy or immunotherapy for intermediate- to high-risk groups. Nevertheless, these are linked with elevated rates of recurrence and disease progression. Within five years of diagnosis, up to 80% of patients with non-muscle invasive bladder cancer (NMIBC) experience relapse, and 30% of these patients exhibit disease progression. According to Zhu et al.[26], 30% of the patient population exhibits muscle invasive bladder cancer (MIBC). According to Cathomas et al.[6], in cases where the disease remains confined to a specific area, potential treatment options for these patient groups may involve radical cystectomy with urinary diversion, radical radiotherapy, and neoadjuvant and/or adjuvant chemotherapy or immunotherapy. According to Powles et al.[7], individuals with advanced metastatic UCB may benefit from cisplatin-containing chemotherapy, which has been demonstrated to enhance overall survival (OS). Additionally, initial findings suggest that immunotherapy utilising PD-1/PD-L1 checkpoint inhibitors may be a promising treatment option.

Immune checkpoint inhibitors (ICIs) have been shown to be a successful therapy option for bladder cancer. Resistance does occur, however, and individuals with bladder cancer respond to ICIs differently[27]. This calls for the investigation of novel therapeutic approaches, with NK cell perhaps developing into a future research hotspot. A growing amount of research points to the possibility that therapies that target NK cell might replace current bladder cancer therapy modalities. First off, one of the key processes for the emergence and growth of tumors is the management of the redox state, which is directly tied to the metabolism of disulfides[28]. For instance, certain tumor cells modify the internal redox environment to promote their survival and tolerance[29]. Disulfides may also have potential use in the therapy of tumors, according to certain research. For instance, cisplatin and paclitaxel, two anti-tumor medications, work by interacting with intracellular disulfides to inhibit tumor growth [30, 31]. SLC7A11 causes NK cell resistance and functions as a regulator in cancers and other illnesses, according to a substantial body of research [32, 33]. In order to increase the therapy response rate and survival rate of bladder cancer patients, it may be necessary to balance NK cell with NK cell.

The study at hand is centred on the examination of genetic variations within NK cell-related genes (DRGs) in relation to bladder cancer. The study involved the analysis of 26 genes linked to NK cell, or DRGs, in a sample of 412 specimens by the researchers. The study revealed that a total of 136 specimens displayed mutations in DRGs, indicating a mutation frequency of 33.01%. The research has identified MYH10 as the gene exhibiting the most frequent mutations, with MYH9, FLNB, and ACTB following suit. The presence of copy number variation (CNV) aberrations was a common occurrence in bladder cancer, particularly in relation to DRGs. The findings of the analysis indicate that a considerable proportion of DRGs manifested copy number deletions, whereas nearly 50% of the DRGs displayed a greater prevalence of CNV amplification. It is worth noting that specific genes, namely SLC3A2, ACTN4, OXSM, and TLN1, demonstrated a greater occurrence of copy number amplification, whereas FLNB, NDUFS1, PDLIM1, and NCKAP1 exhibited a higher frequency of deletions. The research also examined the influence of genetic variations on messenger RNA (mRNA) expression in cases of bladder cancer. The analysis of mRNA expression levels was conducted on both benign and malignant tissues. The findings suggest that Copy Number Variations (CNV) could potentially serve as a mediator for fluctuations in the expression of Dorsal Root Ganglia (DRGs). The DRGs that exhibited a greater frequency of amplification demonstrated a noteworthy upsurge in their expression levels in malignant tissues. Notably, genes such as SLC3A2, OXSM, and GYS1 were found to be overexpressed. In contrast, it was observed that specific genes, namely NDUFS1 and NCKAP1, demonstrated a reduction in their expression levels in cancerous tissues as compared to their non-cancerous counterparts. Moreover, the investigation has identified two discrete phenotypes of disulfidoptosis that are linked to diverse gene expression patterns and characteristics of immune cell infiltration. The results of the cluster analysis indicated the presence of two distinct phenotypes, labelled as NK cell clusters A and B. These phenotypes were found to have divergent prognostic outcomes. The patients categorised within NK cell cluster A exhibited a notably superior prognosis in contrast to those allocated to NK cell cluster B. The cluster B of disulfidoptosis was found to be significantly correlated with the presence of diverse immune cells in the tumour microenvironment, suggesting the possibility of immunosuppressive and inflammatory responses. Enrichment analyses were conducted on differentially expressed genes (DEGs) to gain deeper understanding of the biological behaviours linked with NK cell phenotypes. The tumour microenvironment was found to be implicated in NK cell, as evidenced by the significant enhancement of immune-associated biological pathways. The patients were subjected to unsupervised cluster analysis, resulting in the identification of three distinct clusters that exhibited correlation with NK cell phenotypes. The patients belonging to gene cluster A exhibited a notably advantageous prognosis, whereas those belonging to gene cluster B demonstrated a comparatively disadvantageous prognosis. The computation of a NK cell score was predicated on the variance in gene expression levels that are linked to NK cell. A set of ten genes has been identified as prognostic signatures, out of which six have been categorised as high-risk genes and the remaining four as low-risk genes. The expression levels of certain genes were utilised to calculate the risk score, which exhibited noteworthy distinctions between the NK cell clusters A and B, along with the gene clusters A and B. The predictive capacity of the NK cell score for prognosis was demonstrated across multiple cohorts and subsequently confirmed through validation in testing sets. A prognostic tool for bladder cancer was developed in the form of a nomogram that integrates clinical parameters and risk score.

Despite the remarkable ability of our prognostic signature related to NK cell to characterise the immunological landscape of patients and predict their prognosis, there are still gaps that require attention. Future research must identify appropriate solutions to address these gaps. The utilisation of public databases for data analysis may result in biassed predictions when compared to actual situations. Despite the progress made in mitigating this problem, further input from patients with bladder cancer is necessary to validate the model's suitability and the accuracy of immunotherapy prognosis. In order to refine the relevant details of this investigation, additional prospective studies and fundamental research are necessary.

Declarations of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

We would like to acknowledge all of the subjects for their support with data.

Availability of data and materials

The authors will make the raw data supporting the results of this paper accessible without restriction

Funding

Natural Science Foundation of Hubei Province (2022CFB799)

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, Bray F: Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin 2021, 71(3):209-249.
Siegel RL, Miller KD, Jemal A: Cancer statistics, 2020. CA Cancer J Clin 2020, 70(1):7-30.
Nadal R, Bellmunt J: Management of metastatic bladder cancer. Cancer Treat Rev 2019, 76:10-21.
Alifrangis C, McGovern U, Freeman A, Powles T, Linch M: Molecular and histopathology directed therapy for advanced bladder cancer. Nat Rev Urol 2019, 16(8):465-483.
Choi W, Porten S, Kim S, Willis D, Plimack ER, Hoffman-Censits J, Roth B, Cheng T, Tran M, Lee IL et al: Identification of distinct basal and luminal subtypes of muscle-invasive bladder cancer with different sensitivities to frontline chemotherapy. Cancer Cell 2014, 25(2):152-165.
Cathomas R, Lorch A, Bruins HM, Comperat EM, Cowan NC, Efstathiou JA, Fietkau R, Gakis G, Hernandez V, Espinos EL et al: The 2021 Updated European Association of Urology Guidelines on Metastatic Urothelial Carcinoma. Eur Urol 2022, 81(1):95-103.
Powles T, Kockx M, Rodriguez-Vida A, Duran I, Crabb SJ, Van Der Heijden MS, Szabados B, Pous AF, Gravis G, Herranz UA et al: Clinical efficacy and biomarker analysis of neoadjuvant atezolizumab in operable urothelial carcinoma in the ABACUS trial. Nat Med 2019, 25(11):1706-1714.
Koppula P, Zhuang L, Gan B: Cystine transporter SLC7A11/xCT in cancer: NK cell, nutrient dependency, and cancer therapy. Protein Cell 2021, 12(8):599-620.
Liu X, Nie L, Zhang Y, Yan Y, Wang C, Colic M, Olszewski K, Horbath A, Chen X, Lei G et al: Actin cytoskeleton vulnerability to disulfide stress mediates NK cell. NaNK cell Biol 2023, 25(3):404-414.
Hogg PJ: Biological regulation through protein disulfide bond cleavage. Redox Rep 2002, 7(2):71-77.
Daly EB, Wind T, Jiang XM, Sun L, Hogg PJ: Secretion of phosphoglycerate kinase from tumour cells is controlled by oxygen-sensing hydroxylases. Biochim Biophys Acta 2004, 1691(1):17-22.
Wang Y, Jiang Y, Wei D, Singh P, Yu Y, Lee T, Zhang L, Mandl HK, Piotrowski-Daspit AS, Chen X et al: Nanoparticle-mediated convection-enhanced delivery of a DNA intercalator to gliomas circumvents temozolomide resistance. Nat Biomed Eng 2021, 5(9):1048-1058.
Chen C, Shen M, Liao H, Guo Q, Fu H, Yu J, Duan Y: A paclitaxel and microRNA-124 coloaded stepped cleavable nanosystem against triple negative breast cancer. J Nanobiotechnology 2021, 19(1):55.
Johnson WE, Li C, Rabinovic A: Adjusting batch effects in microarray expression data using empirical Bayes methods. Biostatistics 2007, 8(1):118-127.
Lai H, Cheng X, Liu Q, Luo W, Liu M, Zhang M, Miao J, Ji Z, Lin GN, Song W et al: Single-cell RNA sequencing reveals the epithelial cell heterogeneity and invasive subpopulation in human bladder cancer. Int J Cancer 2021, 149(12):2099-2115.
Wilkerson MD, Hayes DN: ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics 2010, 26(12):1572-1573.
Hanzelmann S, Castelo R, Guinney J: GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics 2013, 14:7.
Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK: limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 2015, 43(7):e47.
Yu G, Wang LG, Han Y, He QY: clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 2012, 16(5):284-287.
Geeleher P, Cox N, Huang RS: pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PLoS One 2014, 9(9):e107468.
Mariathasan S, Turley SJ, Nickles D, Castiglioni A, Yuen K, Wang Y, Kadel EE, III, Koeppen H, Astarita JL, Cubas R et al: TGFbeta attenuates tumour response to PD-L1 blockade by contributing to exclusion of NK cells. Nature 2018, 554(7693):544-548.
Robinson MD, McCarthy DJ, Smyth GK: edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 2010, 26(1):139-140.
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-1756.
Zhang H, Meltzer P, Davis S: RCircos: an R package for Circos 2D track plots. BMC Bioinformatics 2013, 14:244.
Richters A, Aben KKH, Kiemeney L: The global burden of urinary bladder cancer: an update. World J Urol 2020, 38(8):1895-1904.
Zhu CZ, Ting HN, Ng KH, Ong TA: A review on the accuracy of bladder cancer detection methods. J Cancer 2019, 10(17):4038-4044.
Hanna KS, Larson S, Nguyen J, Boudreau J, Bulin J, Rolf M: The role of enfortumab vedotin and sacituzumab govitecan in treatment of advanced bladder cancer. Am J Health Syst Pharm 2022, 79(8):629-635.
Min HY, Lee HY: Oncogene-Driven Metabolic Alterations in Cancer. Biomol Ther (Seoul) 2018, 26(1):45-56.
Iyamu EW: The redox state of the glutathione/glutathione disulfide couple mediates intracellular arginase activation in HCT-116 colon cancer cells. Dig Dis Sci 2010, 55(9):2520-2528.
Sideris S, Aoun F, Zanaty M, Martinez NC, Latifyan S, Awada A, Gil T: Efficacy of weekly paclitaxel treatment as a single agent chemotherapy following first-line cisplatin treatment in urothelial bladder cancer. Mol Clin Oncol 2016, 4(6):1063-1067.
Mitin T, Hunt D, Shipley WU, Kaufman DS, Uzzo R, Wu CL, Buyyounouski MK, Sandler H, Zietman AL: Transurethral surgery and twice-daily radiation plus paclitaxel-cisplatin or fluorouracil-cisplatin with selective bladder preservation and adjuvant chemotherapy for patients with muscle invasive bladder cancer (RTOG 0233): a randomised multicentre phase 2 trial. Lancet Oncol 2013, 14(9):863-872.
Jyotsana N, Ta KT, DelGiorno KE: The Role of Cystine/Glutamate Antiporter SLC7A11/xCT in the Pathophysiology of Cancer. Front Oncol 2022, 12:858462.
Mbah NE, Lyssiotis CA: Metabolic regulation of NK cell in the tumor microenvironment. J Biol Chem 2022, 298(3):101617.

Download PDF

Version 1

posted

You are reading this latest preprint version

ScRNA-seq and bulk RNA-seq reveal the characteristics of NK cell and establish a risk signature in bladder cancer

Status:

Version 1

Abstract

Figures

Introduction

Materials and methods

Bladder Cancer Datasets and Preprocessing

Gene Set Variation Analysis

Estimating of Immune Cell Infiltration

Identification of Differentially Expressed Genes Between NK cell Subtypes in Bladder Cancer

Construction of NK cell Score

Prediction of Multiple Therapeutic Sensitivities

Analysis of ten hub NK cell score related genes by scRNA-seq

Statistical Analyses

Results

Genetic Variations in NK cell-Related Genes among bladder cancer

Secondary Clustering Using Differentially Expressed Genes and Identification of Prognostic-Related Subtypes

The NK cell Score and its Predictive Power for Prognosis in Multiple Cohorts

The NK cell Score and its Associations with Tumor Mutation Burden and Genomic Instability

Immune Cell Infiltration Characteristics and Biological Behaviours between the high risk and low risk group

The NK cell Score Had Great Potential for Predicting Immunotherapy Efficacy

Analysis of scRNA-seq Data of ten hub NK cell score related genes in BLCA

Discussion

Declarations

References

Status:

Version 1