Bioinformatics analysis based on immune-autophagy-related lncRNAs combined with immune infiltration in bladder cancer.

BACKGROUND: To construct a prognostic model based on immune-autophagy-related long noncoding RNA (IArlncRNAs), mainly to predict the overall survival rate (OS) of bladder cancer patients and investigate its possible mechanisms.
METHODS: Transcriptome and clinical data were obtained from The Cancer Genome Atlas (TCGA) database. We identified the IArlncRNA by co-expression analysis, differential expression analysis, and Venn analysis. Then, we identified the independent prognostic IArlncRNAs by univariate Cox regression, LASSO regression, and multivariate Cox regression analysis. Moreover, we constructed the prognostic model based on the independent prognostic IArlncRNAs and clinical features. The proportion of 22 immune cell subtypes was analyzed by the CIBERSORT algorithm. Besides, we identified the differential proportion of 22 immune cell subtypes between the high- and low-risk groups. In addition, we identified the correlation between immune-infiltrating cells (screened by univariate Cox regression and multivariate Cox regression analysis) and IArlncRNAs by Pearson correlation analysis. Finally, we estimated the half-maximal inhibitory concentration (IC50) of chemotherapeutic drugs in patients with bladder cancer based on the pRRophetic algorithm.
RESULTS: Four IArlncRNAs were identified as independent prognostic factors, including AL136084.3, AC006270.1, Z84484.1, and AL513218.1. The OS of patients in the high-risk group was significantly worse compared to the low-risk group. The nomogram showed an excellent predictive effect with the C-index of 0.64. The calibration chart showed a good actual vs. predicted probability. B cells naïve, T cells CD8, T cells CD4 memory resting, T cells follicular helper, macrophages M1, dendritic resting and activated cells had higher infiltrations in the low-risk group and lower infiltration of macrophages M2. The fraction of macrophages M2 was positively associated with AL136084.3. The fraction of T cells CD8 was positively associated with Z84484.1. The fraction of M + macrophages M0 was negatively associated with Z84484.1. Further, we identified the differential IC50 of 24 chemotherapeutic drugs between the high- and low-risk groups.
CONCLUSIONS: The prognostic model based on 4 IArlncRNAs showed an excellent predictive effect. Furthermore, we reasonably speculated that IArlncRNAs are directly or indirectly involved in the immune regulation of the tumor microenvironment (TME), as well as autophagy. 2021 Translational Andrology and Urology. All rights reserved.

Introduction

Bladder cancer is one of the most common malignant tumors of the urinary system, with 430 thousand new cases every year (1). About 25% of the patients with myometrial invasive bladder cancer have a poor prognosis and lack effective treatment (2-4). With the progression of cancer, the 5-years survival rate of patients decreases. In patients with advanced-stage, the mortality rate is more than 50% (5). Therefore, the early evaluation of the prognosis of patients with bladder cancer has important clinical significance.

Clinical features, TMN stage, and lymph node status are critical predictive factors for bladder cancer (6,7). Previous studies have shown that bladder cancer patients with poor clinical conditions, high stage grade, and positive lymph node metastasis have lower OS (8,9). However, traditional predictive methods are mainly based on clinical and pathological information while not considering biological heterogeneity. Recently, new prognostic models have been developed based on glycolysis-related mRNAs, immune-related mRNAs, extracellular matrix-related long noncoding RNA (lncRNAs), microenvironment-related lncRNAs, tumor-infiltrating immune cells, alternative splicing events (10-17).

LncRNA is a non-protein-coding transcript consisting of more than 200 nucleotides (18,19). LncRNAs regulate the occurrence and development of tumors by binding to nuclear acids, proteins, and other macromolecules (20-24). As a result, lncRNAs can also predict cancer patients' outcomes (25,26). Although some molecular biomarkers have been used as independent prognostic factors by traditional predictive models of bladder cancer, there are some limitations, such as small sample size, different detection platforms, single prognostic variates, and similar (6,27). Therefore, it is an urgent clinical task to determine reliable prognostic biomarkers.

The TCGA database works on applying high-throughput genome analysis techniques to better understand the occurrence and progression of carcinoma, thereby improving its diagnosis, treatment, and prevention. Autophagy is a regulatory mechanism that removes unnecessary or dysfunctional cellular components and circulates metabolic substrates. In response to pressure signals in the tumor microenvironment (TME), autophagy pathways in tumor cells and immune cells are altered, thus affecting tumor progression, immunity, and treatment differently. Autophagy is an important mechanism in cancer biology, immunology, and therapy. Thus, in this study, we used a rigorous computational framework to mine immune-autophagy-related long noncoding RNA (IArlncRNAs) and clinical data of bladder cancer downloaded from the TCGA database, as well as to construct a prognostic model IArlncRNAs in order to predict the OS of bladder cancer patients and identify the potential mechanisms. We present the following article in accordance with the TRIPOD reporting checklist (available at https://dx.doi.org/10.21037/tau-21-560).

Methods

Download and organize transcriptional group data and clinical data

The TCGA database is a project supervised by the National Cancer Institute and National Human Genome Research Institute. In this study, we used a TCGA database (https://portal.gdc.cancer.gov/) to download transcriptome (FPKM) and clinical data of patients with bladder cancer (containing primary and metastatic bladder cancer). We organized the transcriptome data and converted the ENSG number to symbol ID. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Identification of IArlncRNAs associated with prognosis of bladder cancer patients

The transcriptome data were divided into mRNAs and lncRNAs. IMMPORT database (https://www.immport.org) was searched for immune-related genes, as well as Human Autophagy Database (http://www.autophagy.lu/) for autophagy genes. We then identified the expression of immune-related and autophagy-related genes from the mRNA dataset by taking intersection with the two target gene sets. We also identified the immune-related lncRNAs and autophagy-related lncRNAs by the co-expression networks by using the R package “limma (version 3.46.0)” (|correlation coefficient| >0.4, and P<0.001) (28). To identify the differential expression of immune-related lncRNAs (DEir-lncRNAs) and autophagy-related (DEar-lncRNAs) lncRNAs, we performed the differential expression analysis by R package limma. Then, IArlncRNAs were identified by taking an intersection between the DEir-lncRNAs and DEar-lncRNAs.

Construction of the prognostic model

Patients who were followed-up for less than 30 days from the clinical data set were excluded. Then, we combined the expression of IArlncRNAs and survival time utilizing the intersection. According to 1:1 ratio, the total samples were randomly divided into the training group and the test group. IArlncRNAs associated with survival were screened in the training group by univariate Cox regression and LASSO regression analysis. The prognostic model was constructed based on the IArlncRNAs above by multivariate Cox regression model, and their coefficients were calculated respectively (we used the R packages “survival” and “survminer”). Then, we calculated the risk score based on the risk coefficient and expression of IArlncRNAs. The formula for calculating the risk score is listed below: Risk score = βgene 1 × exprgene 1 + βgene 2 × exprgene 2 +...+ βgene n × exprgene n.

The patients were divided into the high- and low-risk groups according to the median risk score, after which the heatmaps of IArlncRNAs associated with prognosis in the training group and test group were drawn, respectively. To identify the survival differences between the high- and low-risk groups, we drew the Kaplan-Meier survival curves, survival status distributed plots, and risk score distributed plots. Finally, we verified the predictive value of the risk score model by drawing the ROC curve of the training group and test group (we used the R package “timeROC”).

Establishment of IArlncRNA-mRNA co-expression network

To investigate the association between the IArlncRNAs and their target mRNA, the IArlncRNA-mRNA network was constructed. Pearson correlation coefficients were used to analyze the difference between IArlncRNAs and mRNAs, while the absolute threshold coefficient value >0.4 and P value <0.001 were set as significant. The IArlncRNAs-mRNAs network was visualized by the Cytoscape software (version 3.7.2, http://www.cytoscape.org/).

Independent prognostic analysis and construction of nomogram

We constructed a nomogram by incorporating the related clinical factors (age, gender, stage grade) and risk score obtained from the previous multivariate Cox regression analysis (we used the R package “rms”, “foreign”, and “survival”). We drew the calibration curves and multivariate ROC curves to evaluate the model’s discrimination and accuracy (we used the R package “timeROC”).

CIBERSORT estimation and independent prognostic analysis

To evaluate the proportion of 22 immune-infiltrating cells in each sample, we used a wildly proposed computational algorithm “Cell type Identification By Estimating Relative Subsets Of RNA Transcripts (CIBERSORT)” in R software. The samples with P value <0.05 were considered for further analysis. According to the risk score based on the risk coefficient and expression of IArlncRNAs, we divided the samples into high- and low-risk groups. Then, the Wilcoxon test was used to identify the proportion of immune cell subtypes between the high- and low-risk groups. We identified immune-infiltrating cells associated with survival by univariate Cox regression and multivariate Cox regression (we used the R packages “survival” and “survminer”). To identify the correlation between immune-infiltrating cells and IArlncRNAs, we calculated the Pearson correlation coefficients.

Prediction of the sensitivity of chemotherapeutic agents

To evaluate the differential sensitivity of chemotherapeutic agents between the high- and low-risk groups, pRRophetic algorithm was applied to predict the drugs IC50 based on the Cancer Cell Line Encyclopedia (CCLE) (29,30). R package “pRRophetic” was used to construct the ridge regression model based on the data from the CCLE for the inserted drugs. Besides, we drew boxplots of differential IC50 of interested drugs between the high- and low-risk groups with R software (Wilcoxon signed-rank test).

Statistical analysis

We used R software (version 3.46.0) to perform all statistical analyses, including R packages limma, survival, survminer, timeROC, rms, foreign, and pRRophetic. The prognostic IArlncRNAs were screened by univariate Cox regression, LASSO regression, and multivariate Cox regression analysis. The immune-infiltrating cells associated with survival by univariate Cox regression and multivariate Cox regression. The correlation was determined by Pearson correlation analysis. OS was calculated by the Kaplan-Meier method and evaluated by the log-rank test. The Wilcoxon test was used to identify the proportion of immune cell subtypes between the high- and low-risk groups. Two-tailed P<0.05 was considered statistically significant.

Results

Characteristics of data

Transcriptome data and clinical data on bladder cancer were obtained from the TCGA database. There were 430 specimens in transcriptome data, including 411 tumor specimens and 19 normal specimens. The clinical data of 409 patients included 303 males, 106 females, and 158 deaths. Clinical data showed that the patients’ age ranged from 34 to 90, and their survival and follow-up time ranged from 13 to 5,050 days. Stage grade included 2 stage I, 130 stage II, 139 stage III, and 138 stage IV. These data are shown in .

Table 1

The characteristics of clinical data of bladder cancer in TCGA

ID	Coef	HR	HR.95L	HR.95H	P value
AL136084.3	0.33794746	1.40206684	1.07468	1.82918768	0.01274558
AC006270.1	0.14305056	1.15378814	1.09137872	1.21976638	4.61E-07
Z84484.1	−0.9701779	0.37901561	0.1964704	0.7311678	0.00380431
AL513218.1	−0.524559	0.59181632	0.37154548	0.94267479	0.02720875

Screening of IArlncRNAs

Firstly, we divided the transcriptome data downloaded from the TCGA database into mRNA dataset and lncRNA dataset. Then, we identified the immune-related genes and screened out immune-related lncRNAs by co-expression analysis (correlation coefficient >0.4, and P<0.001), as well as autophagy-related genes and autophagy-related lncRNAs (correlation coefficient >0.4, and P<0.001). Consequently, differential expression analysis of immune-related lncRNAs and autophagy-related lncRNAs between tumor and normal samples was performed (). We obtained 203 DEir-lncRNAs and 121 DEar-lncRNAs. Subsequently, Venn analysis between the DEir-lncRNAs and Dear-lncRNAs was performed, generating 106 IArlncRNAs (). For further analysis, patients with follow-up time <30 days were excluded. The other samples were randomly divided into two groups, including the training group (n=197) and test group (n=196). We further screened out the IArlncRNAs associated with bladder cancer patients’ prognosis in the training group by univariate Cox regression analysis and LASSO regression analysis. We determined a total of 9 IArlncRNAs by univariate Cox regression analysis. The LASSO regression analysis was performed to avoid overfitting of the model. Six IArlncRNAs were more closely related to prognosis (). Subsequently, the multivariate Cox regression analysis was performed by using the significant prognostic factors from the LASSO regression analysis. A total of 4 IArlncRNAs as independent predictors were found, including AL136084.3, AC006270.1, Z84484.1, and AL513218.1 (, ). A correlation coefficient >0 suggested that the expression of lncRNAs were risk factors for the progression of bladder cancer, while the correlation coefficient <0 was interpreted as a protective factor.

Figure 1

Identification of immune-autophagy-related lncRNAs. (A,B) Differential expression analysis of immune-related lncRNAs in the TCGA database; (C,D) differential expression analysis of autophagy-related lncRNAs in the TCGA database. (E) The Venn analysis between DEir-lncRNAs and DEar-lncRNAs. TCGA, The Cancer Genome Atlas; DEir-lncRNAs, differential expression of immune-related lncRNAs; DEar-lncRNAs, differential expression of autophagy-related lncRNAs.

Figure 2

Establishment of a prognostic model of bladder cancer based on IArlncRNAs. (A,B) Based on the univariate Cox regression analysis, the LASSO regression screened 6 IArlncRNAs more closely related to prognosis; (C) the multivariate Cox regression analysis in the training group showed that 4 IArlncRNAs were independent prognostic factors; (D) the overall survival rate of the low-risk group was better than the high-risk group in the training group (P<0.001); (E) the low-risk group’s overall survival rate was better than the high-risk group in the test group (P<0.001); (F) the heatmap showed the differential expression of 4 IArlncRNAs associated with prognosis between the high- and low-risk groups in the training group; (G) the heatmap showed the differential expression of 4 IArlncRNAs associated with prognosis between the high- and low-risk groups in the test group; (H) the distribution of the risk score for each bladder cancer patient in the training group; (I) the distribution of the risk score for each bladder cancer patient in the test group; (J) the survival status of each bladder cancer patient with different risk scores in the training group; (K) the survival status of each bladder cancer patient with different risk scores in the test group; (L) the AUC of 1-year OS was 0.745 (95% CI, 0.660–0.838) in the training group; (M) the AUC of 1-year OS was 0.652 (95% CI, 0.543–0.756) in the test group. IArlncRNAs, immune-autophagy-related lncRNAs; AUC, Area Under Curve; OS, Overall Survival.

Table 2

The IArlncRNAs as independent prognostic factors

Characteristics	Case, n
Age (years)	68.10 (34–90)
Gender
Male	304
Female	106
Tumor stage
T1	3
T2	120
T3	194
T4	59
TX	1
N stage
N0	237
N1	47
N2	76
N3	8
NX	41
M stage
M0	194
M1	11
MX	202
Unknow	2
AJCC stage
Stage I	2
Stage II	130
Stage III	139
Stage IV	136
Unknow	2

A risk score of IArlncRNAs

According to the correlation coefficient and expression of IArlncRNAs, we calculated the risk score for each patient by the formula. The higher the risk score, the worse the prognosis. The patients were divided into the high- and low-risk groups by the median risk score. The results showed that the OS of patients in the high-risk group was significantly worse compared to the low-risk group in the training group () and in the test group (). The heatmap of 4 IArlncRNAs associated with bladder cancer patients’ prognosis showed differential expression of IArlncRNAs between the high- and low-risk groups in the training group () and in the test group (). The risk score of each patient in the training and test groups are shown in . The scatterplot showed the survival status of each bladder cancer patient with a different risk score in the training group () and in the test group ().

To evaluate the sensitivity and specificity of the risk score in predicting the survival of bladder cancer patients, we drew a time-dependent ROC curve. The AUCs of 1-year OS were 0.745 (95% CI, 0.660–0.838) in the training group () and 0.652 (95% CI, 0.543–0.756) in the test group ().

Establishment of the lncRNA-mRNA co-expression network

Subsequently, we identified the potential functions of 4 IArlncRNAs in bladder cancer by establishing the lncRNA-mRNA co-expression network and visualizing it with Cytoscape (). The lncRNA-mRNA co-expression network contained 102 pairs lncRNA-mRNA with | coefficient value | >0.4 and P value <0.05. Among them, 102 mRNAs were significantly associated with 4 IArlncRNAs in the immunity and autophagy signatures. Besides, we drew a Sankey diagram to describe the risk/protective relationship between the 4 IArlncRNAs and 102 mRNAs ().

Figure 3

The construction of IArlncRNA-mRNA co-expression network. (A) The IArlncRNA-mRNA co-expression network contains 102 pairs lncRNA-mRNA with | coefficient value | >0.4 and P value <0.05; (B) the Sankey diagram shows the connection degree between the 4 IArlncRNAs and102 mRNAs (risk/protective). IArlncRNAs, immune-autophagy-related lncRNAs.

Evaluation of 4 IArlncRNAs risk score as independent factors in bladder cancer patients and the construction of the nomogram

Univariate Cox regression and multivariate Cox regression were used to identify whether the age, gender, stage, and risk score were independent prognostic factors of bladder cancer patients. The 95% CI hazard ratio of the risk score in the training group by the univariate Cox regression was 1.038 (1.020–1.056, P<0.001, ), while the 95% CI hazard ratio of risk score by the multivariate Cox regression was 1.036 (1.018–1.054, P<0.001, ). The 95%CI hazard ratio of the risk score in the test group by the univariate Cox regression was 1.092 (0.979–1.217, P=0.113, ), while the 95% CI hazard ratio of risk score by the multivariate Cox regression was 1.144 (1.012–1.293, P=0.032, ). We found that the risk score of 4 IArlncRNAs was an independent factor for bladder cancer patients.

Figure 4

To evaluate the accuracy of IArlncRNAs prognostic signature and clinicopathological features in the prognostic model. (A,B) The forest maps of univariate Cox regression and multivariate Cox regression analysis in the training group; (C,D) the forest maps of univariate Cox regression and multivariate Cox regression analysis in the test group. IArlncRNAs, immune-autophagy-related lncRNAs.

Our nomogram incorporates age, gender, stage, and risk score to predict 1-, 3-, and 5-year overall survival in patients with bladder cancer (). The higher the score, the worse the prognosis. Our results revealed that the risk score had the greatest impact on the prediction of overall survival. In the training group, the C index was 0.64. Also, the calibration curve and ROC curve were drawn (). The prediction was consistent with the actual observation. Similar results were found in the test group ().

Figure 5

The nomogram for the prognostic model. (A) The nomograms of the prognostic model were constructed by including the age, gender, stage, and risk score in the training group; (B) the calibration curves for 1-, 3- and 5-years OS of the nomograms in the training group; (C) the multivariate ROC curves for 1-, 3- and 5-years OS in the training group. The AUC of 1-year OS was 0.746 in the training group. (D) The nomograms of the prognostic model were constructed by including the age, gender, stage, and risk score in the test group; (E) the calibration curves for 1-, 3- and 5-years OS of the nomograms in the test group; (F) the multivariate ROC curves for 1-, 3- and 5-years OS in the test group. The AUC of 1-year OS was 0.654 in the test group. ROC, receiver operating characteristic; AUC, area under curve; OS, overall survival.

Next, we made principal components analysis (PCA) according to IArlncRNAs and total genes’ expression profiles to further understand the difference between the high- and low-risk groups. We found that the risk score can better divide the total samples into two sections, while the cumulative proportion of PC1, PC2, and PC3 were 34.39%, 61.39%, and 83.11% ().

Figure 6

The high- and low-risk patients with bladder cancer showed different immune states. The PCA showed that the risk score could better divide bladder cancer patients into two sections. (A) The PCA between the high- and low-risk groups based on the whole gene expression; (B) the PCA between the high- and low-risk groups based on the 4 IArlncRNAs expressions. PCA, Principal Component Analysis; IArlncRNAs, immune-autophagy-related lncRNAs.

Analysis of the composition and difference of infiltrating immune cells

We used the CIBERSORT algorithm to calculate the constituent proportion of 22 immune cell subtypes in the total samples (). Interestingly, B cells naïve, T cells CD8, T cells CD4 memory resting, T cells follicular helper, macrophages M1, dendritic cells resting, and dendritic cells activated had higher infiltrations in the low-risk groups, while the macrophages M2 infiltration was lower (Wilcoxon test, all P<0.05) (). Subsequently, the immune cell subtypes related to survival were screened by univariate Cox and multivariate Cox regression analysis, while M2 type macrophages and neutrophils were independent prognostic factors ().

Figure 7

The relationship between the IArlncRNAs and infiltrating immune cells. (A) The proportion of 22 immune cell subtypes in bladder cancer patients; (B) the differential proportion of infiltrating immune cells between the high- and low-risk groups; blue and red bar represent the low-risk group and high-risk group, respectively. (C) the result of the multivariate Cox regression for infiltrating immune cells related to survival; (D,E) the correlation between the 4 IArlncRNAs and prognostic infiltrating immune cells. IArlncRNAs, immune-autophagy-related lncRNAs.

We performed co-expression analysis between the 4 prognostic IArlncRNAs and proportions of infiltrating immune cells using the Pearson correlation coefficient (). The fraction of Macrophages M2 was positively associated with AL136084.3 (R=0.38, P=1.4e-05). The fraction of T cells CD8 was positively associated with Z84484.1 (R=0.42, P=9.1e-07). The fraction of macrophages M0 was negatively associated with Z84484.1 (R=−0.23, P=0.011) ().

Differential response of chemotherapy upon high- and low-risk patients

To investigate the differential IC50 of 138 chemotherapeutic drugs between the high- and low-risk groups, pRRophetic algorithm was applied to data from the CCLE. We identified 24 chemotherapeutic drugs containing a variety of anticancer mechanisms (). Besides, we described their targets, whether they have an effect on autophagy of cancer cells, and their references. Further, the boxplots for the differential IC50 of 24 chemotherapeutic drugs between the high- and low-risk groups were drawn (, P<0.05).

Table 3

The characteristics of 24 chemotherapeutic drugs

Drug	Target	Effects on autophagy	References (PMID)
Veliparib (ABT-888)	PARP1 and PARP2 inhibitor	Yes	25975349
Rucaparib (AG-014699)	PARP inhibitor	Yes	31714594
Axitinib	VEGFR-1, VEGFR-2, VEGFR-3, PDGFR-β, and c-KIT	No report	NA
Saracatinib (AZD0530)	Src inhibitor	Yes	32686156
Selumetinib (AZD6244)	MEK inhibitor	Yes	27448918
AZD7762	Chk1 and Chk2 inhibitor	No report	NA
AZD8055	mTOR inhibitor	Yes	31598395
Afatinib (BIBW2992)	EGFR/ErbB inhibitor	Yes	28676644
BMS-754807	IGF-1R/InsR	No report	NA
BX-795	PDK1 and TBK1 inhibitor	Yes	26607461
Erlotinib	EGFR inhibitor	Yes	32190291
Pictilisib (GDC-0941)	PI3K inhibitor	Yes	33048163
Gefitinib	EGFR inhibitor	Yes	33315519
Imatinib	Tyrosine kinase inhibitor	Yes	33066449
JNJ-26854165 (Serdemetan)	HDM2 inhibitor	No report	NA
KU-55933	ATM Kinase Inhibitor	Yes	29921885
Mitomycin C	DNA synthesis inhibitor	Yes	29669822
NU7441	DNA-PK inhibitor	No report	NA
Linsitinib (OSI-906)	IGF-1R inhibitor	No report	NA
Parthenolide	NF-κB inhibitor	Yes	31322140
Pazopanib	VEGFR1, VEGFR2, VEGFR3, PDGFR, FGFR, c-Kit, and c-Fms/CSF1R inhibitor	Yes	23887605
Mirdametinib (PD0325901)	MEK inhibitor	Yes	32268164
PHA-665752	c-Met inhibitor	No report	NA
Roscovitine (Seliciclib)	CDK inhibitor	Yes	30899038

NA, not available.

Figure 8

Differential response of chemotherapy upon high- and low-risk patients. (A-X) A total of 24 chemotherapeutic drugs with differential IC50 between high- and low-risk groups.

Discussion

Bladder cancer is a common malignant tumor of the urinary system with a high recurrence rate and different progression rates. At present, the most common treatments for bladder cancer include surgery, radiotherapy, and chemotherapy. The progressive development of new treatment strategies has significantly improved the survival rate of bladder cancer patients. However, the risk of mortality in patients with advanced bladder cancer is still high, and there is a lack of an excellent prognostic prediction model. At present, postoperative follow-up monitoring of the recurrence of bladder cancer patients mainly depends on cystoscopy. Yet, this invasive method limits its application in large-scale cancer surveillance. In addition, urine cytology is poorly used in predicting low-grade bladder cancer. Identifying biomarkers of bladder cancer may improve the ability to screen and diagnose bladder cancer and judge its prognosis. Therefore, there is an urgent requirement to explore the novel and effective prognostic indicators of bladder cancer.

Bioinformatics technology is a powerful high-throughput tool for screening molecular biomarkers and prognostic indicators. LncRNA is a non-protein-coding transcript with a length of more than 200 nucleotides and different functions in different physiological and pathological conditions. LncRNAs have tissue-specific expression patterns and a vital role in various diseases, including cancer. For example, Zhan et al. (31) found that SOX2OT is positively correlated with histological grade, TNM stage, and poor prognosis and may serve as a potential prognostic factor in bladder cancer. Moreover, Zheng et al. (32) suggested lncRNA PTENP1 as a promising biomarker for detecting bladder cancer; lncRNA PTENP1 has a suppressor role in cancer progression and promotes prognosis. Besides, Zhan et al. (33) also showed an interesting result with a great reference value that PCAT-1 acts as an independent prognostic factor for recurrence-free survival of bladder cancer and presents a potential value for clinical diagnosis and prognosis.

In the present study, we identified 4 IArlncRNAs with immunity and autophagy signatures as independent factors that contributed to the survival of bladder cancer patients, including AL136084.3, AC006270.1, Z84484.1, and AL513218.1. According to the risk score, the patients were divided into high- and low-risk groups and analyzed the two groups’ immune infiltration. We found that the proportion of B cells naïve, T cells CD8, T cells CD4 memory resting, T cells follicular helper, macrophages M1, dendritic resting and activated cells, and macrophages M2 were different between the two groups, which may contribute to the regulation of immunity. Endogenous tumor-associated antigens (TAA), recognized and presented to T lymphocytes by dendritic cells (DCs) with the major histocompatibility complex (MHC), leads to the adaptive anti-cancer immune responses (34,35). Classic MHC includes MHC I, MHC II, and MHC III; when MHC I presents self-proteins and viral proteins to CD8+ T lymphocytes, MHC II presents extracellular antigens to CD4+ T lymphocytes. In many types of tumor cells, MHC reduction is often observed, which enables tumor cell immune escape, resulting in tumor progression and tolerance to immunotherapy (36).

Autophagy has a dual role in antigen presentation, the key process of anti-cancer immune responses. Activation of autophagy promotes antigen presentation to CD8+T lymphocytes mediated by DCs, which then stimulates cytotoxic responses (37,38). Besides, autophagy can also disrupt the process of antigen presentation. As an autophagy cargo receptor gene, NBR1 targets and degrades the MHC I, thus disrupting its antigen presentation ability to CD8+ T cells, which in turn can be reversed by autophagy inhibitor (39). Meanwhile, mucin domain protein-4 (TIM-4) binds to AMPK-a1 and activates the autophagy, and degrades the TAA, thus disrupting the antigen presentation and leading to the reduction of CD8+ T cells (40).

TME, which consists of the tumor and other cells (stromal cells, extracellular matrix (ECM), and immune cells, such as natural killer (NK) cells, tumor-associated macrophages (TAMs), myeloid-derived suppressor cells (MDSCs), and T and B lymphocytes), is essential for the progression, maintenance, metastasis, and response to therapy of different types of cancers (41). Previous studies have shown that the immunosuppressive compartments of TME, including the autophagy activity of MDSCs, M2 macrophage polarization, and T cells regulatory (Tregs) recruitment, are associated with immune evasion (13,42,43). In our study, CD8 T cells and DCs had higher, while the M2 macrophage had lower infiltrations in the low-risk groups. Therefore, we reasonably speculated that IArlncRNAs are directly or indirectly involved in the immune regulation of TME, as well as autophagy.

This study has a few limitations. First, the data obtained from public databases were limited, and the clinical data were not comprehensive, which may lead to potential errors or deviations. Second, in this study, we did not verify the stability of IArlncRNAs prognostic characteristics in another independent data set and further experiments. Third, all the data used in this study were obtained from Western countries; therefore, careful consideration should be taken when applying our conclusions to patients from Asian countries.

Conclusions

In this study, we constructed a prognostic model with excellent predictive ability based on 4 IArlncRNAs in patients with bladder cancer. IArlncRNAs are critical factors, directly or indirectly involved in the immune regulation of TME and autophagy.

The article’s supplementary files as