A Novel Assessment Model Based on Molecular Subtypes of Hypoxia-Related LncRNAs for Prognosis of Bladder Cancer

Hypoxia is a common feature in various tumors that regulates aggressiveness. Previous studies have demonstrated that some dysregulated long non-coding RNAs (lncRNAs) are correlated with tumor progression, including bladder cancer (BCa). However, the prognostic effect of hypoxia-related lncRNAs (HRLs) and their clinical relevance, as well as their regulatory effect on the tumor immune microenvironment, are largely unknown in BCa. A co-expression analysis between hypoxia genes and lncRNA expression, which was downloaded from the TCGA database, was performed to identify HRLs. Univariate Cox regression analysis was performed to select the most desirable lncRNAs for molecular subtype, and further LASSO analysis was performed to develop a prognostic model. This molecular subtype based on four HRLs (AC104653, AL136084, AL139393, and LINC00892) showed good performance in the tumor microenvironment and tumor mutation burden. The prognostic risk model suggested better performance in predicting BCa patients’ prognosis and obtained a close correlation with clinicopathologic features. Furthermore, four of five first-line clinical chemotherapies showed different sensitivities to this model, and nine immune checkpoints showed different expression in the molecular subtypes or the risk model. In conclusion, this study indicates that this molecular subtype and risk model based on HRLs may be useful in improving the prognostic prediction of BCa patients with different clinical situations and may help to find a useful target for tumor therapy.

Introduction

More than 570,000 patients were diagnosed with bladder cancer (BCa) in the past year, making BCa one of the top 10 cancers with high morbidity and mortality in the world, which is a common malignant tumor of the urinary system with a poor prognosis (Sung et al., 2021). In current clinical practice, pathological biopsy with cystoscopy is considered to be the most reliable method for detecting this highly heterogeneous cancer (Cumberbatch and Noon, 2019). In all BCa cases, more than 75% of cases are non-muscle invasive (NMIBC), and the primary and traditional treatment option mainly includes surgery (transurethral resection of bladder tumors, TURBT) with or without chemotherapy (intravesical therapy, IVT) (Bryan et al., 2021). However, almost a quarter of NMIBC patients ultimately develop muscle-invasive BCa (MIBC) because of its high recurrence rate and the occurrence of chemical resistance (Babjuk et al., 2019). In recent years, although immunotherapies using immune checkpoint inhibitors have gradually been used as a treatment, there are no biomarkers to assess the effectiveness of immunotherapy in BCa. Therefore, novel and accurate therapeutic evaluations for BCa treatment are urgently needed.

Cellular hypoxia is a common feature in various physiological or pathophysiological states in the tumor microenvironment (TME) and plays critical roles in cancer initiation, progression, and treatment resistance in various cancers, including BCa (Reiher et al., 2001; Wan et al., 2014). With the development of gene-sequencing technology, many hypoxia-related signatures have been established in multiple tumors to predict patient prognosis and treatment efficacy, such as immunotherapy and chemotherapy response (Shao et al., 2021; Sun et al., 2021a). In BCa, multiple hypoxia-related gene signatures were developed and validated to predict prognosis and TME immune characteristics (Yang et al., 2017; Liu et al., 2021). In recent years, with the successful application of immunotherapy in tumors, an increasing number of people have begun to pay attention to the clinical value of immune checkpoints, such as anti-PD-L1/PD-1 therapy. It is worth noting that the effect of immunotherapy is closely related to the TME (Hu et al., 2021; Liu et al., 2021). Liu et al. (2021) developed and validated a hypoxia risk score that was related to the expression of several immune checkpoints, such as PD-L1, PD-1, CTLA-4, and LAG-3, by analyzing the BCa TME.

Long non-coding RNAs are untranslated RNAs of more than 200 nucleotides in length, and over 15,000 lncRNAs are encoded by the human genome (Zhang et al., 2017). Dysregulated lncRNA expression is correlated with cancer progression and has been reported in different cancer types, including BCa (Liang et al., 2018). Several hypoxia-related lncRNAs (HRLs) have been certified to play crucial roles in tumor regulation or processes (Yuan et al., 2019). For instance, hypoxic cells remodel the tumor microenvironment to facilitate tumor growth and development by secreting oncogenic lncRNA UCA1 enriched exosomes (Xue et al., 2017). LncRNA HOXA-AS3 could act as a competing endogenous RNA (ceRNA) of miR-455-5p to regulate Notch1 and regulate chemotherapeutic drug sensitivity in BCa cells (Chen et al., 2020).

However, only a few studies have reported the role of several HRLs in BCa without its clinical prognostic value and its impact in the TME. Therefore, understanding HRL differential expression in BCa patients may help to identify biomarkers that can act as useful therapeutic targets and predict patient prognosis.

The present study intends to classify the expression of HRLs in BCa patients and develop a model with bioinformatics methods to investigate the influence of the TME, and predict patient prognosis, chemotherapy, and immunotherapy response.

Materials and Methods

Data Source

The expression of RNAseq (including mRNAs and lncRNAs) and clinical information (including age, gender, tumor grade, disease stage, TNM stage, and prognosis) were obtained from The Cancer Genome Atlas database (TCGA¹) by using the keywords “TCGA-BLCA AND HTSeq-FPKM AND Gene Expression Quantification” on April 4, 2021. Patients with incomplete information were excluded.

The hypoxia genes were searched from the Molecular Signatures Database v7.4 of GSEA². The gene was excluded if gene expression was not detected in more than half of the BCa patient samples.

Correlation Analysis

To identify the HRLs, Pearson correlation analysis was performed between different expression levels of hypoxia genes (log₂FC of > 1 and FDR of < 0.05) and all lncRNA expression data in samples to identify the HRLs (| correlation coefficient| ≥ 0.4 and p < 0.001). The heatmap of correlation analysis was performed using the OmicStudio tools at https://www.omicstudio.cn/tool.

Cell Lines and RT-qPCR Analysis

A normal human bladder epithelial cell line (SV-HUC1) and BCa cell lines (T24, TCCSUP, and UMUC-3) were obtained from the American Type Culture Collection (ATCC; Manassas, VA, United States). The cells were cultured in RPMI 1640 medium (Hyclone), supplemented with 10% fetal bovine serum (Gibco) and 1% penicillin/streptomycin, in a humidified incubator (Thermo) at 37°C with 5% CO₂. Total RNAs were extracted from the collection cells using Trizol Kit (Invitrogen, Carlsbad, CA, United States) according to the manufacturer’s instructions. The first strand (cDNA) was reverse-transcribed and used as the template for qPCR analysis as previously described (Wang et al., 2020). The expression levels of lncRNAs were measured using the SYBR Green qPCR Kit (Takara, Japan). The gene name and primer sequences were listed in Supplementary Table 1. Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) was used as an endogenous control.

Consensus Clustering

Univariate Cox regression analysis was used to extract the HRLs via the R package “survival” (p-value ≤ 0.05). Consensus clustering was performed to explore the potential molecular subtype between the BCa patients using the R package “ConsensusClusterPlus.” The count of cluster (k) was set from 1 to 9, and the best optimal k value was selected for further investigation.

Gene Set Enrichment Analysis (GSEA)

Gene set enrichment analysis (version 4.1.0) was performed to identify the differences in the set of genes expressed between the two clusters in the enrichment of the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway. The number of permutations was performed 1,000 times for each analysis, and the pathways with FDR < 0.05 were considered statistically significant.

Estimation of Tumor Immune Microenvironment

An analytical tool, CIBERSORT, was used to evaluate the tumor immune infiltration levels of BCa by providing an estimation of the abundances of 22 kinds of cell types based on all gene expression levels. The algorithm was run for 1,000 permutations and BCa samples with an output p < 0.05 were selected for further analysis. The immune score, stromal score, and ESTIMATE score were calculated by the R package “ESTIMATE.”

Estimation of Tumor Mutational Burden (TMB)

Tumor mutational burden is an emerging therapeutic measure of sensitivity to immunotherapy. The R package “maftools” was used to process somatic mutation data including somatic, coding, base replacement, and insert–deletion mutations. The median TMB value was defined as the cutoff value to separate BCa patients into High-TMB and Low-TMB groups.

Development of the Hypoxia-Related LncRNAs Prognostic Model

All BCa samples were randomly equally divided into construction and verification datasets at a ratio of 1:1. A least absolute shrinkage and selection operator (LASSO) Cox hazards model analysis was used to construct the HRL risk model with the optimal prognostic value by the R package “glmnet.” The following formula was used to calculate the risk score of every patient according to the regression coefficients from the model and lncRNAs expression: $Riskscore={\sum }_i^n(coef\times expression)$ (n: the number of lncRNAs; i: the serial number of each lncRNA; coef: the coefficient value from the Cox hazards analysis; expression: HRL expression). Univariate and multivariate analyses were performed to verify the independent prognostic predictors via the R package “survival.” A nomogram was used to predict the 1-, 3-, and 5-year survival ratio based on the clinicopathologic features (such as age, gender, tumor grade, and disease stage) and risk score of this model to assist clinical procedures and improve risk stratification by using the R package “rms.”

Predicting Chemotherapeutic Response

Based on the information retrieved from the Genomics of Drug Sensitivity in Cancer (GDSC) database³, the chemotherapeutic sensitivity for each BCa was calculated. Five common clinical chemotherapeutic drugs were selected to predict the chemotherapeutic response, and the R package “pRRophetic” was used to estimate the chemotherapeutic response determined by the half-maximal inhibitory concentration (IC50).

Statistical Analysis

R software (version 4.0.3), and GraphPad Prism (version 6.0) were used to analyze the statistical significance. One-way ANOVA followed by ad hoc Tukey’s multiple comparisons was utilized to compare the expression of lncRNAs in different cell lines. The Wilcoxon signed-rank test was used to analyze continuous variables, whereas the chi-square test was used to analyze categorical data. Kaplan–Meier curves were plotted, and a log-rank test was used to calculate the significant survival difference. The sensitivity and specificity of prognostic prediction using the HRL risk score were analyzed by Receiver Operating Characteristic (ROC) analysis. The area under the ROC curve (AUC) served as an indicator to validate prognostic accuracy. A p < 0.05 was considered statistically significant.

Results

Identification of Hypoxia-Related LncRNAs

A total of 16,808 lncRNAs and 20,454 mRNAs from 414 BCa tumor samples and 19 non-tumor samples were downloaded from the TCGA database. Of the 198 hypoxia genes selected from the GSEA database, 65 differentially expressed hypoxia genes were identified in the TCGA-BLCA cohort (Supplementary Figures 1A,B). Among them, the expression matrixes of 15 hub hypoxia-related genes were extracted. In total, by co-expression analyses, 67 HRLs were confirmed (Supplementary Figures 1C,D). Combined with patient prognostic information downloaded from TCGA datasets, univariate Cox regression was then implemented to screen HRLs from these 67 genes (Figures 1A,B), and four HRLs (AC104653, AL136084, AL139393, and LINC00892) were significantly related to BCa patient prognosis and significantly differentially expressed in tumor and non-tumor specimens (Figures 1C,D). Correspondingly, RT-qPCR analysis of basic experiments also demonstrated that the expression level of four HRLs was lower in tumor cell lines (T24, TCCSUP, and UMUC-3) than in a normal bladder cell line (SV-HUC1) (Figure 1E).

FIGURE 1

Figure 1. Identification of a molecular subtype of HRLs. (A) Study flow chart. (B) The heatmap of the correlations between hypoxia-related genes and hypoxia-related lncRNAs. (C) The four hypoxia-related prognostic lncRNAs. (D) The expression of four lncRNAs. (E) The expression level of four HRLs in a normal bladder cell line (SV-HUC1) and tumor cell lines (T24, TCCSUP, and UMUC-3), *p < 0.05, **p < 0.01 and ***p < 0.001.

The Molecular Subtype Based on the Hypoxia-Related LncRNAs

To investigate the connections between the expression of these four HRLs and BCa subtypes, consensus clustering analysis was performed with BCa patients. The results showed that when k = 2, patients are divided into two clusters with the highest intragroup correlations and the lowest intergroup correlations (Supplementary Figure 2). A relationship between patient clusters and the gene expression profiles and their clinical features, such as grade, stage, TNM stage, gender, and age, are shown in a heatmap, indicating that there is a significant difference in disease stage between the two clusters (Figure 2A). Additionally, the Kaplan–Meier curve showed that patients in Cluster 1 had a significantly prolonged overall survival (OS) ratio compared with those in Cluster 2 (Figure 2B).

FIGURE 2

Figure 2. The molecular subtype based on the HRLs. (A) The heatmap of the correlations of molecular subtype and clinical features. (B) Kaplan–Meier curve of two molecular subtypes, *p < 0.05, **p < 0.01, and ***p < 0.001.

Gene Set Enrichment Analysis Enrichment Analysis of Molecular Subtypes

To investigate which pathways were enriched in clusters, GSEA was performed to identify the significant KEGG pathways associated with the two clusters in BCa patients. Patients in Cluster 2 were mainly enriched in the MAPK signaling pathway, TGFβ signaling pathway, angiotensin system, ligand receptor interaction, smooth muscle contraction, calcium signaling pathway, leukocyte transendothelial migration, and focal adhesion. However, patients in Cluster 1 were enriched in oxidative phosphorylation (FDR < 0.05) (Supplementary Figure 3).

Immune Function Analysis of the Molecular Subtypes

To explore the difference in immune function in the two clusters, the ESTIMATE score, immune score, and stromal score of BCa samples were calculated. The results showed that patients in Cluster 2 had a significantly higher score than those in Cluster 1 in three kinds of scores (Figures 3A–C). Further investigation of the specific distribution of immune cell fractions in BCa indicated that the Cluster 1 group contained a higher proportion of follicular helper T cells while the Cluster 2 group contained a higher proportion of T cells CD8 and mast cells resting (Figure 3D). Meanwhile, more than 80 immune genes were differentially expressed in these two clusters, indicating that HRL molecular subtypes are associated with immune genes and the tumor microenvironment (Supplementary Figure 4). On the other hand, in recent years, with the advancement of the concept of precision medicine, immune checkpoint inhibitor (ICI) therapy has attracted great attention. Therefore, the different expression levels of ICIs in different subtypes were further explored. Among nine known immune checkpoints (CTLA4, PD-L1, PD-L2, PD-1, SIGLEC15, TIM3, TIGIT, IDO1, and LAG3), there was no difference between normal and tumor samples (Supplementary Figure 5). Remarkably, the overall expression levels of genes including CTLA4, TIM3, and PD-L2 were relatively higher in Cluster 2 while SIGLEC15 was higher in Cluster 1, and the other five genes showed no difference in expression (Figure 3E).

FIGURE 3

Figure 3. Estimation of the function of the immune cells. The boxplot of correlation between molecular subtypes and immune score (A), stromal score (B), and estimate score (C). (D) Relationships between the molecular subtype and the immune cell infiltration in BCa. (E) The expression of nine immune checkpoints in different molecular subtypes, *p < 0.05, **p < 0.01, and ***p < 0.001.

Tumor Mutational Burden Analysis of the Molecular Subtypes

After BCa patients were divided into H-TMB (n = 113) and L-TMB (n = 209) by the median value of the TMB score, the waterfall plots revealed the top mutated genes in the two clusters (Figures 4A,B). Overall, comparing the most frequent somatic mutations between the two clusters, TP53 and TTN were similar between the two clusters in terms of their mutation frequencies. However, the mutation frequencies of ARID1A, KMT2D, and RB1 were ranked third, fourth, and fifth in Cluster 2 (33, 30, and 30%, respectively) while the mutation frequencies were lower than 30% in Cluster 1. The Kaplan–Meier results indicated that H-TMB can prolong patient survival compared with the L-TMB (Figure 4C). Further investigation in cluster subtypes showed that patients in Cluster 2 had a lower TMB score than those in Cluster 1 (Figure 4D). Interestingly, in a study on whether patients in different clusters with different TMB scores would have different prognoses, the Kaplan–Meier analysis demonstrated that patients in Cluster 1 with H-TMB had a significantly better prognosis than patients in other subtypes (Figure 4E).

FIGURE 4

Figure 4. Estimation of the tumor mutation burden (TMB). The waterfall plot of mutation genes in cluster 1 (A) and cluster 2 (B). (C) The boxplot of the TMB score in the different molecular subtypes. (D) Kaplan–Meier curve of H-TMB and L-TMB. (E) Kaplan–Meier curve of different TMB scores and molecular subtypes.

Development of a Prognostic Risk Model

Bladder cancer patients who matched their clinical features and survival information were selected from all BCa samples. Based on four tumor classification lncRNAs, three optimal lncRNAs (AC104653, AL136084, and LINC00892) were selected to construct a model due to the coefficient caused by the LASSO Cox regression analysis (Figures 5A,B). The risk score was calculated according to the following formula: risk score = (−0.005 × AC104653 expression) + (0.431 × AL136084 expression) + (−1.361 × LINC00892 expression). According to this model, all BCa patients were divided into two risk groups with the median score as the cutoff. Principal component analysis (PCA) was used to further verify the risk model, and the results confirmed complete separation between different risk sample sets in the PCA distribution 3D plot (Figures 5C–F). In addition, the Kaplan–Meier analysis of patient prognosis showed that patients in the low-risk group had a higher OS ratio than those in the high-risk group (Supplementary Figure 6).

FIGURE 5

Figure 5. Construction of the prognostic risk model. (A,B) LASSO regression of the minimum criteria calculation. PCA of risk score based on the expression profiles of all genes (C), hypoxia genes (D), all HRLs (E), and HRL prognostic risk-related genes (F).

Validation of the Risk Model

To reveal the accuracy of the risk score, BCa samples were randomly divided into two cohorts for construction and verification. In both sets of cohorts, the distribution of the risk scores and vital statuses of patients, sorted by the risk score, was the same (Figures 6A,B). The Kaplan–Meier curve results also suggested that the OS ratio of high-risk patients was significantly lower than that of the low-risk group (Figures 6C,D). Moreover, time-dependent ROC analyses also found that this model was an effective way to predict BCa patient prognosis (Figures 6E,F). Univariate and multivariate analyses were used to further investigate BCa patients’ independent prognostic factors and indicated that the risk score was an independent prognostic factor in the construction cohort, considering patient age, gender, and clinical stage as in the verification cohort (Figures 6G–J). Based on the ROC analysis results generated by the risk model, the AUC value was 0.68, which was higher than that of clinicopathologic features, including age, gender, grade, and disease stage, indicating that the prognostic model has a better ability to predict the survival of patients (Supplementary Figure 7).

FIGURE 6

Figure 6. The construction and verification of the HRL risk model of BCa patients. Distribution of the three lncRNA-based risk scores and vital statuses of patients in different risk groups in the construction (A) and verification cohort (B). Kaplan–Meier survival curves of OS of different risk patients in construction (C) and verification cohort (D). ROC curve of the prognostic model in the construction (E) and verification cohort (F). (G,H) The univariate and multivariate analyses of patients in the construction cohort. (I,J) The univariate and multivariate analyses of patients in the verification cohort.

Stratification Analysis of the Risk Score Clinical Utility

To further explore the clinical value of this risk score with different clinical features, BCa patients were separated into different groups based on age, gender, grade, and disease stage. Considering the different stratified analysis results, the BCa patients with high-risk scores predict significantly poorer OS than those with low-risk scores aged ≤ 60 or > 60. Likewise, the results were confirmed in females or males, with high tumor stage (T1–T3), with high disease stage (stage II–IV), with high tumor grade, and without lymph node or distant metastasis (N0 and M0) (Figure 7A). Of note, the results to identify whether clinicopathological features were associated with the risk score revealed a significantly increasing trend in age, T stage, and N stage (Supplementary Figure 8). Referring to all the above findings, based on the risk model score and clinicopathological features, the nomogram was established as a clinically applicable quantitative tool to predict the OS of BCa patients (Figure 7B). In addition, considering the clinical treatment value, the high-risk patients showed more sensitivity to gemcitabine, doxorubicin, methotrexate, and vinblastine than the low-risk group, while there was no significant difference in sensitivity to cisplatin (Figure 7C).

FIGURE 7

Figure 7. The clinical utility of the risk score. (A) Kaplan–Meier survival curves of OS of high and low risk in different clinical stratification. (B) Prognostic survival is predicted in the nomographic chart. (C) Different chemotherapeutic responses in different risk BCa patients.

Associations Between the Risk Model and Molecular Subtype and Immune Function

A Sankey diagram was established to illustrate the relationship of molecular subtype and risk score and patient prognosis. The results showed that more than half of the patients in Cluster 2 could be classified as high-risk and predict poor survival (Figure 8A). Moreover, the boxplot results also supported that BCa patients’ risk score in Cluster 2 was significantly higher than that of patients in Cluster 1 (Figure 8B). For the correlation between immune score and risk score, the boxplot showed that a high immune score was negatively correlated with a risk score (Figure 8C). The heatmap analysis also confirmed this point (Figure 8D). All the above findings suggested that patients in Cluster 2 were more likely to obtain high-risk scores and poor OS and that there was a negative correlation between immune score and risk score. To further investigate the relationship between immune cell infiltration and the risk score, the results revealed that the risk score was negatively related to the functions of CD4 memory activated T cells, follicular helper T cells, regulatory T cells, and CD8 T cells, while positively related with functions of mast cells activated, macrophages M0, and macrophages M2 (Figure 8E). In addition, the correlation between the risk score and the immune checkpoint was also analyzed and demonstrated that the risk score was negatively related to the IDO1, LAG3, PD-1, PD-L1, CTLA4, and TIGIT gene expression (Figure 8F).

FIGURE 8

Figure 8. The relationship between the risk model and molecular subtype and immune function. (A) The Sankey diagram of the relationship between molecular subtype, risk score, and prognosis. (B) The boxplot of the relationship between risk score and molecular subtypes. (C) The boxplot of the relationship between risk score and immune score groups. (D) The heatmap of the correlations of risk score with molecular subtype and clinical features. (E) Correlations between risk score and immune infiltrating level of seven types of immune cells. (F) Comparison of the expression of selected immune checkpoint in high- and low-risk groups.

Discussion

As one of the most common tumors of the urinary system, the main treatment for BCa in current clinical practice includes surgical removal of the tumor, such as TURBT, and supplementation with chemotherapy, such as cisplatin. However, with tumor recurrence and drug resistance, BCa patients’ prognosis has not improved within the past several years (Antoni et al., 2017). Recently, immunotherapy as a novel treatment method for multiple cancers has received increasing attention. A clinical trial demonstrated that patients who received atezolizumab had a better prognosis and manageable side effects than patients who received chemotherapy (van der Heijden et al., 2021). However, due to the lack of effective inclusion criteria for the use of immunotherapy, immunotherapy has not become the first choice of clinical first-line treatment.

One of the main reasons why tumors are different from normal tissues is that tumor cells will continue to adapt to the bad environment where they are. Hypoxia is a poor environment that weakens the function of the tumor, while malignant tumor cells are often able to compensate for the process and drive the occurrence of more malignant disease behaviors later (Chae et al., 2016). Therefore, hypoxia as a factor that alters the tumor characteristic microenvironment can accelerate tumor metastasis and enhance drug resistance (Zandberg et al., 2021). Numerous studies have confirmed that there is a hypoxic microenvironment in most solid tumors, including BCa, and the expression of hypoxia-related genes will affect the TME, the efficacy of therapeutic drugs, and patient prognosis (Yang et al., 2017; Liu et al., 2021). However, some previous studies have verified that some HRLs will also affect tumor progression, and the HRL risk model can also predict the patient prognosis well and has an effect on clinical therapy (Shao et al., 2021; Zhang et al., 2021).

In the present study, BCa patients were separated into two molecular subtypes, and the immune microenvironment was evaluated. Univariate Cox regression was used to screen HRLs for subtyping. Four HRLs’ (AC104653, AL136084, AL139393, and LINC00892) low expression levels were used to construct the subtype model. GSEA showed hypoxia genes involved in the MAPK signaling pathway, and the TGFβ signaling pathway was more enriched in the Cluster 2 group than in Cluster 1. TMB analysis showed that different clusters had significantly different mutated genes, such as RB1, KDM6A, and ARID1A, as well as prognosis, which needs further research to recognize its gene functions. Further investigation of the immune microenvironment showed patients in Cluster 2, which had higher immune scores, had higher proportions of CD8 T cells and mast resting cell phenotypes. Moreover, more than 80 differentially expressed immune genes also indicated that hypoxia caused by tumor cell proliferation and malformation of the tumor vascular systems can significantly affect the TME. In addition, based on this subtype, three of them (AC104653, AL136084, and LINC00892) were selected to establish a risk score model by LASSO regression, and patients who obtained higher risk scores predicted bad outcomes compared with the low-risk score. Interestingly, high-risk score patients are sensitive to multiple chemotherapy drugs except for cisplatin compared with low-risk patients. Moreover, nine immune checkpoint analyses also showed that four genes (CTLA4, PD-L2, SIGLEC15, and TIM3) based on these two clusters showed significant difference, and six genes (TIGIT, IDO1, LAG3, PD-1, PD-L1, and CTLA4) have been shown significant differences based on this risk score.

Hypoxia-inducible factor-1 alpha (HIF-1α) has been shown to be activated in the hypoxic microenvironment and to regulate the TGFβ signaling pathway by targeting TGFβ1 (Lazzara et al., 2020). Hypoxia can induce DNA hydroxymethylation and increase the expression of TNFα in the MAPK signaling pathway (Wu et al., 2015). These results are consistent with the results of the GSEA in this study.

Tumor mutational burden is a significant biological marker to indicate the state of tumor mutation and has been described as a powerful tool for identifying candidates associated with tumor immune regulation mechanisms (Klebanov et al., 2019; Cai et al., 2021). Remarkably, in the current study, although the survival analysis illustrated that L-TMB has a worse prognosis than H-TMB, patients in Cluster 2 had a worse survival rate whether in L-TMB or H-TMB compared with Cluster 1, suggesting that patients in Cluster 2 may be most in urgent need of targeted therapies or combination therapies.

The functions of immune cell types are influenced by hypoxia, thereby directly or indirectly inducing tumor development and serving as a barrier to immunotherapy. The expression of HIF-1α in the hypoxic microenvironment promotes cancer glycolysis and evasion of immunosurveillance to promote the function of CD8⁺ T cells (Gropper et al., 2017). In addition, CD4⁺ T cells and NK cells are also regulated to affect anti-tumor immunity and checkpoint molecule expression in the hypoxic microenvironment (Westendorf et al., 2017; Ou et al., 2019). Therefore, these studies have demonstrated that different immune cell subpopulations could respond to hypoxia quite differently, which further confirms the results of this study.

The rise of biological big data mining for patient prognosis has promoted the rapid development of precision medicine in recent years (Mollayeva et al., 2019). Previous studies have performed HRL analyses for multiple types of tumors and built prognostic models, including clear cell renal carcinoma (Zhang et al., 2021). However, there is no further research or model establishment in BCa. Currently, a risk score model is established by LASSO regression, and patients in Cluster 2 obtain higher risk scores than those in Cluster 1. The ROC curve also can illustrate that this kind of risk model in predicting prognosis has more sensitivity and specificity than other clinical features, such as age, gender, tumor grade, or TNM stage. Meanwhile, patients with different risk scores show different outcomes in different clinical groups.

Immune checkpoint inhibitors are a potential treatment for cancer therapy that block key molecules and have shown impressive efficacy against cancer, especially revolutionizing the management of metastatic and locally advanced BCa (de Kouchkovsky et al., 2021; Sun et al., 2021b). In the multifaceted immune regulation of BCa, CTLA4 is one of the primary immune checkpoints. CTLA4 is highly increased in the tumor immune microenvironment, particularly on T regulatory cells (Tregs) upon activation (Arce Vargas et al., 2018). Ipilimumab (IgG1 mAb) is an approved ICI drug of the anti-CTLA4 antibody to be used as a treatment of advanced melanoma and has demonstrated a significant and long-lasting effect in approximately 15 to 20% of treated patients (Schadendorf et al., 2015). PD-L2, similar to PD-L1, is another ligand of PD-1, that has a prominent role in maintaining self-tolerance under normal physiological conditions, limiting T cell activation and proliferation in peripheral tissues (Dowell et al., 2021). Previous studies have demonstrated that multiple tumor patients with PD-L2 expression had worse outcomes, even in the absence of PD-L1, including BCa (Gao et al., 2009; Ariafar et al., 2020). SIGLEC15 is an emerging broad-spectrum target for normalization cancer immunotherapy and is complementary to PD-L1 (Wang et al., 2019). SIGLEC15 can inhibit the CD8⁺ T cell proliferation to promote tumor growth, and the clinical trials of its inhibitors are also ongoing (Hu et al., 2021). TIM3 is expressed on several immune cell types, such as CD4⁺ Th1 cells, CD8⁺ T cells, Tregs, NK cells, and monocytes (Monney et al., 2002; Gorman and Colgan, 2014). However, unlike other immune checkpoints, the ligands for TIM3 have been expressed on a variety of cancer cells (Gorman and Colgan, 2014). In addition, the anti-TIM3 antibodies are being developed (Rotte et al., 2018). TIGIT is expressed in both T cells and NK cells and inhibits immune response mediated via triggering CD155 on tumor cells (Anderson et al., 2016; Attalla et al., 2020). LAG3 is expressed on both natural and induced Tregs, T cells, B cells, and DCs, and its receptor ligands include MHC class II molecules expressed on APCs and Galectin-3 (Rotte et al., 2018). IDO1 can inhibit effector T cells and hyper-activating Tregs to induce immune suppression in T cells, while it is not expected to kill tumor cells directly (Jung et al., 2019). However, an early report showed that the combined use of IDO1 inhibitors and PD-1 or PD-L1 inhibitors may enhance the efficacy of checkpoint blockade alone (Jung et al., 2019). In the current study, these nine immune checkpoints had a significant difference in molecular subtype or risk score, suggesting that this molecular subtype and risk score may provide the likelihood of immunotherapy for different BCa patients.

Based on clinical treatment guidelines, cisplatin-based chemotherapy remains the standard therapy in perioperative and first-line metastatic settings (Patel et al., 2020). For patients who have been deemed cisplatin-resistant and harbor tumors with high PD-L1 expression, PD-L1 inhibitor drugs would be used for the first line (Patel et al., 2020). Luckily, according to this study, combining the results of risk scores and molecular subtypes, different patients can choose more sensitive chemotherapy drugs and targeted immunotherapy drugs, single or combined. This is more in line with the concept of individualized treatment in precision medicine.

Although the HRL model has been developed in multiple tumors, this study is the first to develop an HRL model in BCa with some experimental verification. These HRLs were developed based on data from the TCGA database, which contains complete gene expression and clinical and prognosis data of patients with BCa. Meanwhile, the basic experiment in BCa cell lines also verified the expression of HRLs, which is consistent with the result from the TCGA database. However, the current study still has some limitations. First, after searching the GEO database, which contains many gene sequence cohorts, no other available independent lncRNA cohort can be used to validate the usefulness of this model in BCa. Meanwhile, because of the lack of BCa samples, this study cannot verify its effect in clinical practice. In addition, the underlying mechanisms of these four lncRNAs should be investigated by many specific experiments.

Conclusion

The molecular subtype and risk assessment model based on four HRLs (AC104653, AL136084, AL139393, and LINC00892) can improve the prognostic prediction of BCa patients with different clinical situations. Furthermore, the different expression of ICIs and the different sensitivity of five chemo drugs included in this molecular subtype and risk model might be beneficial targets to explore the development of BCa and design personalized individualized therapy strategies.

Data Availability Statement

Publicly available datasets were analyzed in this study. This data can be found here: https://portal.gdc.cancer.gov/, http://www.gsea-msigdb.org/gsea/.

Author Contributions

XJ and XC constructed this study. XC and YZ collected and analyzed the data and drafted the manuscript. FW, XZ, and QF made the figures. XY, JL, and XC revised the manuscript. All authors read and approved the final manuscript.

Funding

This work was supported by the Natural Science Foundation of Zhejiang Province, Zhejiang Province, China (Grant No. LGF18H050001), the Scientific Funding of Health Commission of Zhejiang Province (Grant No. 2019KY420), and the National Natural Science Foundation of China (NSFC Grant No. 81900595).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

Thanks to the technical support from the Central Laboratory, the First Affiliated Hospital, Zhejiang University School of Medicine.

Supplementary Material

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

Supplementary Figure 1 | Identification of hypoxia genes and lncRNAs. The heatmap (A) and a volcano plot (B) of the different expressions of hypoxia genes between normal and BCa samples. (C) The protein-protein interaction (PPI) network of 15 hub hypoxia genes. (D) The mRNA-lncRNA network of hypoxia genes and HRLs.

Supplementary Figure 2 | The heatmap of the consensus clustering of HRLs.

Supplementary Figure 3 | GSEA for the different molecular subtypes of HRLs.

Supplementary Figure 4 | The heatmap of the expression of different molecular subtypes.

Supplementary Figure 5 | The boxplot of the expression of 9 immune checkpoints in normal and BCa samples.

Supplementary Figure 6 | The Kaplan-Meier curve of OS of the risk model.

Supplementary Figure 7 | The Survival-dependent ROC curves of risk model prognostic indicators.

Supplementary Figure 8 | The boxplot of the clinicopathological significance of the risk model.

Footnotes

1. ^ https://portal.gdc.cancer.gov/
2. ^ http://www.gsea-msigdb.org/gsea/
3. ^ https://www.cancerrxgene.org/

References

Anderson, A. C., Joller, N., and Kuchroo, V. K. (2016). Lag-3, Tim-3, and TIGIT: co-inhibitory receptors with specialized functions in immune regulation. Immunity 44, 989–1004. doi: 10.1016/j.immuni.2016.05.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Antoni, S., Ferlay, J., Soerjomataram, I., Znaor, A., Jemal, A., and Bray, F. (2017). Bladder cancer incidence and mortality: a global overview and recent trends. Eur. Urol. 71, 96–108. doi: 10.1016/j.eururo.2016.06.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Arce Vargas, F., Furness, A. J. S., Litchfield, K., Joshi, K., Rosenthal, R., Ghorani, E., et al. (2018). Fc effector function contributes to the activity of human anti-CTLA-4 antibodies. Cancer Cell 33, 649–663.e4.

Google Scholar

Ariafar, A., Ghaedi, M., Rezaeifard, S., Shahriari, S., Zeighami, S., Ghaderi, A., et al. (2020). Clinical relevance and prognostic significance of PD-1/PD-Ls in non-metastatic bladder cancer: a role for PD-L2. Mol. Immunol. 124, 35–41. doi: 10.1016/j.molimm.2020.05.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Attalla, K., Farkas, A. M., Anastos, H., Audenet, F., Galsky, M. D., Bhardwaj, N., et al. (2020). TIM-3 and TIGIT are possible immune checkpoint targets in patients with bladder cancer. Urol. Oncol. (in press). doi: 10.1016/j.urolonc.2020.06.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Babjuk, M., Burger, M., Comperat, E. M., Gontero, P., Mostafid, A. H., Palou, J., et al. (2019). European association of urology guidelines on non-muscle-invasive bladder cancer (TaT1 and carcinoma in situ) – 2019 update. Eur. Urol. 76, 639–657. doi: 10.1016/j.eururo.2019.08.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Bryan, R. T., Liu, W., Pirrie, S. J., Amir, R., Gallagher, J., Hughes, A. I., et al. (2021). Comparing an imaging-guided pathway with the standard pathway for staging muscle-invasive bladder cancer: preliminary data from the bladderpath study. Eur. Urol. 80, 12–15. doi: 10.1016/j.eururo.2021.02.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Cai, Y., Wang, X., Wang, N., Wu, J., Ma, L., Xie, X., et al. (2021). Correlations between tumor mutation burden and immune infiltrates and their prognostic value in pancreatic cancer by bioinformatic analysis. Life Sci. 277:119505. doi: 10.1016/j.lfs.2021.119505

PubMed Abstract | CrossRef Full Text | Google Scholar

Chae, Y. C., Vaira, V., Caino, M. C., Tang, H. Y., Seo, J. H., Kossenkov, A. V., et al. (2016). Mitochondrial Akt regulation of hypoxic tumor reprogramming. Cancer Cell 30, 257–272. doi: 10.1016/j.ccell.2016.07.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, D., Xie, S., Wu, Y., Cui, Y., Cai, Y., Lan, L., et al. (2020). Reduction of bladder cancer chemosensitivity induced by the effect of HOXA-AS3 as a ceRNA for miR-455-5p that upregulates notch1. Front. Oncol. 10:572672. doi: 10.3389/fonc.2020.572672

PubMed Abstract | CrossRef Full Text | Google Scholar

Cumberbatch, M. G. K., and Noon, A. P. (2019). Epidemiology, aetiology and screening of bladder cancer. Transl. Androl. Urol. 8, 5–11. doi: 10.21037/tau.2018.09.11

PubMed Abstract | CrossRef Full Text | Google Scholar

de Kouchkovsky, I., Zhang, L., Philip, E. J., Wright, F., Kim, D. M., Natesan, D., et al. (2021). TERT promoter mutations and other prognostic factors in patients with advanced urothelial carcinoma treated with an immune checkpoint inhibitor. J. Immunother. Cancer 9:e002127. doi: 10.1136/jitc-2020-002127

PubMed Abstract | CrossRef Full Text | Google Scholar

Dowell, A. C., Munford, H., Goel, A., Gordon, N. S., James, N. D., Cheng, K. K., et al. (2021). PD-L2 is constitutively expressed in normal and malignant urothelium. Front. Oncol. 11:626748. doi: 10.3389/fonc.2021.626748

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, Q., Wang, X. Y., Qiu, S. J., Yamato, I., Sho, M., Nakajima, Y., et al. (2009). Overexpression of PD-L1 significantly associates with tumor aggressiveness and postoperative recurrence in human hepatocellular carcinoma. Clin. Cancer Res. 15, 971–979. doi: 10.1158/1078-0432.CCR-08-1608

PubMed Abstract | CrossRef Full Text | Google Scholar

Gorman, J. V., and Colgan, J. D. (2014). Regulation of T cell responses by the receptor molecule Tim-3. Immunol. Res. 59, 56–65. doi: 10.1007/s12026-014-8524-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Gropper, Y., Feferman, T., Shalit, T., Salame, T. M., Porat, Z., and Shakhar, G. (2017). Culturing CTLs under hypoxic conditions enhances their cytolysis and improves their anti-tumor function. Cell Rep. 20, 2547–2555. doi: 10.1016/j.celrep.2017.08.071

PubMed Abstract | CrossRef Full Text | Google Scholar

Hu, J., Yu, A., Othmane, B., Qiu, D., Li, H., Li, C., et al. (2021). Siglec15 shapes a non-inflamed tumor microenvironment and predicts the molecular subtype in bladder cancer. Theranostics 11, 3089–3108. doi: 10.7150/thno.53649

PubMed Abstract | CrossRef Full Text | Google Scholar

Jung, K. H., LoRusso, P., Burris, H., Gordon, M., Bang, Y. J., Hellmann, M. D., et al. (2019). Phase I study of the indoleamine 2,3-dioxygenase 1 (IDO1) inhibitor navoximod (GDC-0919) administered with PD-L1 inhibitor (Atezolizumab) in advanced solid tumors. Clin. Cancer Res. 25, 3220–3228. doi: 10.1158/1078-0432.CCR-18-2740

PubMed Abstract | CrossRef Full Text | Google Scholar

Klebanov, N., Artomov, M., Goggins, W. B., Daly, E., Daly, M. J., and Tsao, H. (2019). Burden of unique and low prevalence somatic mutations correlates with cancer survival. Sci. Rep. 9:4848. doi: 10.1038/s41598-019-41015-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Lazzara, F., Trotta, M. C., Platania, C. B. M., D’Amico, M., Petrillo, F., Galdiero, M., et al. (2020). Stabilization of HIF-1alpha in human retinal endothelial cells modulates expression of miRNAs and proangiogenic growth factors. Front. Pharmacol. 11:1063. doi: 10.3389/fphar.2020.01063

PubMed Abstract | CrossRef Full Text | Google Scholar

Liang, H., Zhao, X., Wang, C., Sun, J., Chen, Y., Wang, G., et al. (2018). Systematic analyses reveal long non-coding RNA (PTAF)-mediated promotion of EMT and invasion-metastasis in serous ovarian cancer. Mol. Cancer 17:96. doi: 10.1186/s12943-018-0844-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Z., Tang, Q., Qi, T., Othmane, B., Yang, Z., Chen, J., et al. (2021). A robust hypoxia risk score predicts the clinical outcomes and tumor microenvironment immune characters in bladder cancer. Front. Immunol. 12:725223. doi: 10.3389/fimmu.2021.725223

PubMed Abstract | CrossRef Full Text | Google Scholar

Mollayeva, T., Sutton, M., Chan, V., Colantonio, A., Jana, S., and Escobar, M. (2019). Data mining to understand health status preceding traumatic brain injury. Sci. Rep. 9:5574. doi: 10.1038/s41598-019-41916-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Monney, L., Sabatos, C. A., Gaglia, J. L., Ryu, A., Waldner, H., Chernova, T., et al. (2002). Th1-specific cell surface protein Tim-3 regulates macrophage activation and severity of an autoimmune disease. Nature 415, 536–541. doi: 10.1038/415536a

PubMed Abstract | CrossRef Full Text | Google Scholar

Ou, Z. L., Luo, Z., Wei, W., Liang, S., Gao, T. L., and Lu, Y. B. (2019). Hypoxia-induced shedding of MICA and HIF1A-mediated immune escape of pancreatic cancer cells from NK cells: role of circ_0000977/miR-153 axis. RNA Biol. 16, 1592–1603. doi: 10.1080/15476286.2019.1649585

PubMed Abstract | CrossRef Full Text | Google Scholar

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: 10.3322/caac.21631

PubMed Abstract | CrossRef Full Text | Google Scholar

Reiher, F. K., Ivanovich, M., Huang, H., Smith, N. D., Bouck, N. P., and Campbell, S. C. (2001). The role of hypoxia and p53 in the regulation of angiogenesis in bladder cancer. J. Urol. 165(6 Pt 1), 2075–2081. doi: 10.1016/S0022-5347(05)66297-0

CrossRef Full Text | Google Scholar

Rotte, A., Jin, J. Y., and Lemaire, V. (2018). Mechanistic overview of immune checkpoints to support the rational design of their combinations in cancer immunotherapy. Ann. Oncol. 29, 71–83. doi: 10.1093/annonc/mdx686

PubMed Abstract | CrossRef Full Text | Google Scholar

Schadendorf, D., Hodi, F. S., Robert, C., Weber, J. S., Margolin, K., Hamid, O., et al. (2015). Pooled analysis of long-term survival data from phase II and phase III trials of ipilimumab in unresectable or metastatic melanoma. J. Clin. Oncol. 33, 1889–1894. doi: 10.1200/JCO.2014.56.2736

PubMed Abstract | CrossRef Full Text | Google Scholar

Shao, J., Zhang, B., Kuai, L., and Li, Q. (2021). Integrated analysis of hypoxia-associated lncRNA signature to predict prognosis and immune microenvironment of lung adenocarcinoma patients. Bioengineered 12, 6186–6200. doi: 10.1080/21655979.2021.1973874

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, X., Luo, H., Han, C., Zhang, Y., and Yan, C. (2021a). Identification of a hypoxia-related molecular classification and hypoxic tumor microenvironment signature for predicting the prognosis of patients with triple-negative breast cancer. Front. Oncol. 11:700062.

Google Scholar

Sun, X., Zhou, Z., Zhang, Y., Wang, J., Zhao, X., Jin, L., et al. (2021b). Identification and validation of a hypoxia-related prognostic and immune microenvironment signature in bladder cancer. Cancer Cell Int. 21:251.

Google Scholar

Sung, H., Ferlay, J., Siegel, R. L., Laversanne, M., Soerjomataram, I., Jemal, A., et al. (2021). Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J. Clin. 71, 209–249. doi: 10.3322/caac.21660

PubMed Abstract | CrossRef Full Text | Google Scholar

van der Heijden, M. S., Loriot, Y., Duran, I., Ravaud, A., Retz, M., Vogelzang, N. J., et al. (2021). Atezolizumab versus chemotherapy in patients with platinum-treated locally advanced or metastatic urothelial carcinoma: a long-term overall survival and safety update from the phase 3 IMvigor211 clinical trial. Eur. Urol. 80, 7–11. doi: 10.1016/j.eururo.2021.03.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Wan, G., Xie, W., Liu, Z., Xu, W., Lao, Y., Huang, N., et al. (2014). Hypoxia-induced MIR155 is a potent autophagy inducer by targeting multiple players in the MTOR pathway. Autophagy 10, 70–79. doi: 10.4161/auto.26534

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, F., Fan, M., Cai, Y., Zhou, X., Tai, S., Yu, Y., et al. (2020). Circular RNA circRIMS1 acts as a sponge of miR-433-3p to promote bladder cancer progression by regulating CCAR1 expression. Mol. Ther. Nucleic Acids 22, 815–831. doi: 10.1016/j.omtn.2020.10.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J., Sun, J., Liu, L. N., Flies, D. B., Nie, X., Toki, M., et al. (2019). Siglec-15 as an immune suppressor and potential target for normalization cancer immunotherapy. Nat. Med. 25, 656–666. doi: 10.1038/s41591-019-0374-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Westendorf, A. M., Skibbe, K., Adamczyk, A., Buer, J., Geffers, R., Hansen, W., et al. (2017). Hypoxia enhances immunosuppression by inhibiting CD4+ effector T cell function and promoting treg activity. Cell. Physiol. Biochem. 41, 1271–1284. doi: 10.1159/000464429

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, M. Z., Chen, S. F., Nieh, S., Benner, C., Ger, L. P., Jan, C. I., et al. (2015). Hypoxia drives breast tumor malignancy through a TET-TNFalpha-p38-MAPK signaling axis. Cancer Res. 75, 3912–3924. doi: 10.1158/0008-5472.CAN-14-3208

PubMed Abstract | CrossRef Full Text | Google Scholar

Xue, M., Chen, W., Xiang, A., Wang, R., Chen, H., Pan, J., et al. (2017). Hypoxic exosomes facilitate bladder tumor growth and development through transferring long non-coding RNA-UCA1. Mol. Cancer 16:143. doi: 10.1186/s12943-017-0714-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, L., Taylor, J., Eustace, A., Irlam, J. J., Denley, H., Hoskin, P. J., et al. (2017). A gene signature for selecting benefit from hypoxia modification of radiotherapy for high-risk bladder cancer patients. Clin. Cancer Res. 23, 4761–4768. doi: 10.1158/1078-0432.CCR-17-0038

PubMed Abstract | CrossRef Full Text | Google Scholar

Yuan, S., Xiang, Y., Wang, G., Zhou, M., Meng, G., Liu, Q., et al. (2019). Hypoxia-sensitive LINC01436 is regulated by E2F6 and acts as an oncogene by targeting miR-30a-3p in non-small cell lung cancer. Mol. Oncol. 13, 840–856. doi: 10.1002/1878-0261.12437

PubMed Abstract | CrossRef Full Text | Google Scholar

Zandberg, D. P., Menk, A. V., Velez, M., Normolle, D., DePeaux, K., Liu, A., et al. (2021). Tumor hypoxia is associated with resistance to PD-1 blockade in squamous cell carcinoma of the head and neck. J. Immunother. Cancer 9:e002088. doi: 10.1136/jitc-2020-002088

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, H., Qin, C., Liu, H. W., Guo, X., and Gan, H. (2021). An effective hypoxia-related long non-coding RNAs assessment model for prognosis of clear cell renal carcinoma. Front. Oncol. 11:616722.

Google Scholar

Zhang, P., He, D., Xu, Y., Hou, J., Pan, B. F., Wang, Y., et al. (2017). Genome-wide identification and differential analysis of translational initiation. Nat. Commun. 8:1749. doi: 10.1038/s41467-017-01981-8

PubMed Abstract | CrossRef Full Text | Google Scholar