The RNA sequencing (RNA-seq) expression profile data, somatic mutation data, and survival data regarding pan-cancer (33 species) were taken from the database of UCSC Xena (https://xenabrowser.net/). The format of RNA-seq data was changed from Fragments Per Kilobase Million (FPKM) to the format of Transcripts per million (TPM), and then we did a log2 conversion. Among them, analysis and processing of the downloaded somatic mutation data were done by mutect. Finally, the copy number variations (CNV) data processed by the gistic algorithm was also provided by the UCSC Xena database (http://xena.ucsc.edu/), while the methylation data was taken from the LinkedOmics database (http://linkedomics.org).

The BLCA GEO queue was retrieved from the GEO database (https://www.ncbi.nlm.nih.gov/geo/), which has extensive survival data, such as GSE31684, GSE48075, GSE13507, GSE32894, GSE48277, and GSE69795. The BLCA samples were kept.

Afterward, we also downloaded three cohorts linked with immunotherapy, GSE78220 (melanoma), GSE135222 (NSCLC), and GSE91061 (melanoma). Following the knowledge sharing 3.0 License Agreement, the complete expression data and comprehensive clinical data of the IMvigor210 queue (BLCA immunotherapy-related queue) came from http://research-pub.Gene.com/imvigor210corebiologies/.

The analysis difference between subtypes was done using the limma package (19), and differentially expressed genes were chosen through the | log2 (Fold Change) | >1 and False Discovery Rate (FDR) <0.05.

We enriched the differentially expressed genes among subtypes and then carried out an analysis using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) through the WebGestaltR package (version 0.4.2) (20), and the chosen gene set was “c2 cp. kegg. v7.0. symbols. Gmt”, which had the KEGG channel. The GSEA input file consisted of the expression profile data. The threshold values of enriched pathways were p<0.05 and FDR < 0.25. Likewise, the GO function enrichment analysis of differentially expressed genes was done using the R software package WebGestaltR (threshold value was set as P < 0.05).

The R-package survminer (https://cran.r-project.org/package=survminer) was employed to get the best cutoff of genes from various datasets, and the samples were sorted into high and low expression groups following the best cutoff, and afterward, we drew the KM curve. We randomly collected the BLCA cancer samples from the TCGA dataset using the ratio train:test = 7:3. Univariate analysis was done in the training data set. The R software package glmnet (21) was utilized for establishing the Lasso expression model (COX). Based on the model created in this study, the most suitable genes were chosen using single factor Cox regression, and we obtained 21 genes when the value of minimum lambda=0.04090851, which were employed for multivariate analysis. StepAIC approach was utilized for reducing gene number. The stepwise regression used the Akaike information criterion (AIC) (22), which took into account the model’s statistical fit and the number of parameters used for fitting. The stepAIC procedure in the MASS package started from the most complex model and successively eliminated a variable to lower the AIC. The smaller the value, the better the model, which suggested that the model had a sufficient fit with fewer parameters. 11 genes were obtained at the end according to this procedure.

Through a comprehensive study of hundreds of various tumor expression profiles, the Tumor Immune Dysfunction and Exclusion (TIDE) (http://TIDE.dfci.harvard.edu/) (23) analysis can uncover biomarkers that predict the therapeutic response of immune checkpoint inhibitors/medicines. The TIDE score obtained from TIDE analysis can be used to determine the sensitivity of immunological checkpoints.

TMB is a quantifiable immune-response biomarker that reflects the number of mutations in tumor cells. TMB scores were calculated using a Perl script and corrected by dividing by the total length of exons.

Base on genes from previous research (24) and ssGSEA analysis was used to analyze these genes to define T cell inflamed score.

The difference in clinicopathological features among the three subtypes was investigated using the Chi-square test. The expression levels of three subtypes were determined using ANOVA. The difference in the two groups was investigated using the T-test. For correlation analysis, the Pearson correlation coefficient was used. R (version 4.0.2) was used for statistical analyses. Statistical significance was defined as a P-value of < 0.05.

We discovered four types of genes in the literature including MHC, chemokine, immune-stimulator, and receptor (25). The Spearman correlation between these genes and the ISCA1 was measured in pan-cancer. The link between these four gene types with ISCA1 was varied in various types of cancer. It was mostly positive in uveal melanoma (UVM), BLCA, kidney papillary cell carcinoma (KIRP), etc. while thyroid carcinoma (THCA), testicular germ cell tumors (TGCT), etc. were mostly linked negatively ( Figure 1A ). Moreover, the link between CTLA4, PDCD1, CD86, CD274, and ISCA1 in various cancer types was measured. The outcomes revealed that these four genes were substantially positively linked with ISCA1 in BLCA ( Figures 1B–E ). Furthermore, ssGSEA was employed for evaluating the scores of 28 immune cell scores in various cancer types, and then their link with ISCA1 was measured. The outcomes of this analysis suggested that there was a considerable positive link between the expression of ISCA1 and 20 immune scores in BLCA ( Figure 1F ).

Moreover, the ISCA1 gene expression in pan-cancer was observed ( Figure 2 ). The outcomes indicated that: among the 24 tumors with para-cancerous samples, the expression of the ISCA1 gene in 15 cancer species was considerably varied in comparison to that in para-cancerous samples. Among them, the expression of the ISCA1 gene was lowered in tumor samples of 11 cancer species, including breast invasive carcinoma (BRCA), BLCA, cervical squamous cell carcinoma (CSCC), glioblastoma multiform (GBM), kidney chromophobe (KICH), endocervical adenocarcinoma (CESC), kidney renal clear cell carcinoma (KIRC), KIRP, THCA, endometrial carcinoma (UCEC), lung adenocarcinoma (LUAD), rectum adenocarcinoma (READ), etc. Moreover, the expression of cholangiocarcinoma (CHOL), head and neck squamous cell carcinoma (HNSC), liver cancer (LIHC), and stomach cancer (STAD) were enhanced in tumor samples of four cancer species.

Based on the above analysis of pan-cancer, the ISCA1 gene had a considerable positive link with four types of genes in BLCA: immune-stimulator, MHC, chemokine, and receptor. In BLCA, the ISCA1 gene was positively linked with CTLA4, PDCD1, CD86, CD274, and immune score. Using the difference analysis, it was discovered that the expression of the ISCA1 gene in BLCA tumor samples was decreased in comparison to that in adjacent samples. Furthermore, the survival analysis highlighted that samples were divided into the high ISCA1 expression group and low ISCA1 expression group., so we were focused on the role of ISCA1 in BLCA.

In herein, we used the surv_cutpoint function of the SurvMiner package to find the best cutoff for grouping. BLCA samples were categorized into two groups ( Figure 3A ) as per the best cutoff of ISCA1 expression value. The group with high ISCA1 expression showed a poor prognosis, suggesting that ISCA1 is a risk factor in BLCA. Afterward, we mapped 10 genes with the highest mutation frequencies in the high and low expression groups. The outcomes revealed that TTN, synb1, TP53, RB1, kmt2d, arid1a, and other genes in the low expression group had lower mutation frequencies ( Figure 3B ). However, no major variation was observed in the TMB of high and low expression groups of ISCA1 ( Figure 3C ). We observed that there was a major difference in the expression of ISCA1 with CNV amplification and deletion and the normal copy number, and the expression of ISCA1 with CNV amplification was greatly enhanced, while that with CNV deletion was greatly reduced (p < 0.0001, Figure 3D ). Meanwhile, no major link was seen between the methylation degree and the expression of ISCA1 ( Figure 3E ). Also, using limma analysis, 1672 genes, included 1505 upregulated genes and 167 downregulated genes, were screened between high ISCA1 group and low ISCA1 group ( Figure S1A ). KEGG pathways enrichments analysis showed that 1505 genes were enriched in 10 KEGG pathways, such as, PI3K-Akt signaling pathway, cell cycle ( Figure S1B ), while 167 genes were not enriched into the KEGG pathway.

In the ISCA1 expression group, the differential expression of chemokine, immune-stimulator, MHC, and receptor genes was investigated ( Figures 4A–E ). The expression was higher in the high expression group, and most of the four different types of genes had substantial variations.

The ISCA1 high expression group had a high immune score, and the distribution differences of 28 immune cell scores in the ISCA1 group were examined, revealing that 21 had major variations ( Figure 5A ). The immune infiltration of BLCA was next investigated, and a link was observed between ISCA1 expression and immune infiltration score. The marker genes of five cell types were studied: CD8 T cells, dendritic cells, macrophages, NK cells, and Th1 cells. In the ISCA1 high expression group, the majority of the genes were highly expressed ( Figure 5B ). The link between ISCA1 and immunological checkpoints was also measured. Based on the outcomes, ISCA1 and these immunological checkpoints had a substantial positive link ( Figure 5C ).

The link between ISCA1 expression value and pan-cancer T cell inflamed score was measured, and the outcomes indicated a major positive link ( Figure 6A ). Moreover, the link between ISCA1 and different immune properties (immune checkpoint, expression of immunomodulator and TIIC effector genes, and characteristics linked with immunotherapy) in subgroups with varied immune checkpoint blockade (ICB) responses were analyzed ( Figure 6B ). The outcomes suggested that ISCA1 had a major positive link with them.

By comparing the scores of BLCA tumors and immune-related pathways, it was discovered that there were major variations in related immune pathways in the high ISCA1 group and the low ISCA1 group of BLCA tumors, the correlation between ISCA1 expression and Neuroendocrine differentiation pathway is positive in low ISCA1 group ( Figure 6C ). For example, the Neuroendocrine_ differentiation pathway score was higher in the high ISCA1 group, whereas the score was lower in the low ISCA1 group. ARID1A, RB1, ERBB2, FANCC, and other genes that could be linked to radiotherapy and chemotherapy were compared. It was observed that the mutation frequencies in the high and low ISCA1 groups were different ( Figure 6D ). ARID1A, RB1, ERBB2, ERCC2, and FANCC mutation frequencies were greater in high ISCA1 groups, for instance, missense mutation of ERCC2 was 5% in the high ISCA1 group and 4% in the low ISCA1 group. In high and low ISCA1 expression groups, the differences in three categories (EGFR network, immune inhibit oncogenic pathways, and radiotherapy predicted pathways) were compared ( Figure 6E ). The high ISCA1 group was mostly positively correlated in the EGFR network in the EGFR_ligands pathway, while the low ISCA1 group was mostly negatively linked.

A total of 575 up-regulated genes and 100 down-regulated genes were obtained by grouping the up-regulated and down-regulated genes of the BLCA sample species ISCA1, StromalScore, and ImmuneScore, respectively ( Figures 7A, B ). The GO and KEGG function enrichment of differential genes was then analyzed using WebGestaltR. Genes were discovered to be closely linked to cancer and immune pathways such as myeloid leukocyte migration, leukocyte migration, angiogenesis, Th1 and Th2 cell differentiation, and so on ( Figures 7C–F ).

Using the string website, PPI found and analyzed a total of 675 differential genes. Following that, Cytoscape was used to visualize the data and the MCODE plug-in was utilized to identify significant clusters. There were three gene clusters with more than ten genes each (Mcode1, Mcode2, and Mcode3 respectively). The genes MRC1, CXCL11, CCL3, CCL4, CSF1, and FN1 were all found in Mcode1 ( Figure 8A ). Then, to determine their functions, WebGestaltR was utilized to do a GO and KEGG function enrichment analysis ( Figures 8B–E ). The findings revealed that the Mcode1 module was linked to immunological pathways such as the Toll-like receptor signaling pathway and the interaction between cytokine and cytokine receptors.

After the above analysis, 675 we identified the differential genes linked with immunity, and then 172 genes linked with prognosis were provided by univariate analysis (p < 0.05). Afterward, Lasso was employed for selecting the most suited gene for developing the IRS model. Based on the minimum lambda = 0.04090851, we obtained 21 genes ( Figure 9A ). These genes were used for the multivariate analysis. To further decrease gene number, the stepAIC approach was employed. Finally, we got 11 genes ( Figure 9B ), and the risk coefficients of linked genes were obtained. The risk scores of each sample in the training and validation datasets were measured, and the best cutoff score was used to categorize them into high and low-risk groups, with their KM curves and ROC curves demonstrated separately. In the training set, the AUC value for the 1-year survival rate was 0.81, the AUC value for the 3-year survival rate was 0.75, and the AUC value for the 5-year survival rate was 0.77, whereas in the test set, the AUC value for the 1-year survival rate was 0.75, the AUC value for the 3-year survival rate was 0.72, and the AUC value for the 5-year survival rate was 0.64. ( Figures 9C, D ). A greater survival rate (p < 0.0001) was observed in the low-risk group in both the training and validation sets. In addition, all TCGA datasets, GSE13507 datasets, and GSE32894 datasets were used to validate our IRS model ( Figure 9E–G ). The aforesaid datasets’ BLCA tumor samples could also be sorted into high- and low-risk groups with differing prognoses (p < 0.0001), with the low-risk group having a greater rate of survival.

The immunotherapy datasets IMvigor210, GSE91061, GSE78220, and GSE135222 were chosen to predict, evaluate, and compare the efficacy scores of immune therapy. Our approach was used to calculate IRS in these data, and TIDE was utilized to evaluate the effect of immunotherapy, after which the predictive effect of IRS and TIDE on treatment response was evaluated. The immunotherapy samples were separated into high and low-scoring groups following the best IRS and TIDE cut-off scores. Our IRS score was higher than the TIDE score, according to the results ( Figure 10 ).