Construction and analysis of a ceRNA network and patterns of immune infiltration in bladder cancer.

BACKGROUND: Bladder cancer (BC) is the ninth most common malignant tumor, accounting for an estimate of 549,000 new BC cases and 200,000 BC-related deaths worldwide in 2018. The prognosis of BC has not substantially improved despite significant advances in the diagnosis and treatment of the disease.
METHODS: The RNA sequencing (RNA-seq) data and clinical information of BC patients were downloaded from The Cancer Genome Atlas (TCGA) database. The Cell-type Identification By Estimating Relative Subsets Of RNA Transcripts (CIBERSORT) algorithm was used to assess immune infiltration. The survival analyses were performed using the selected components of a ceRNA network and selected immune cell types by least absolute shrinkage and selection operator (LASSO) Cox regression to calculate the risk score. The accuracy of prognosis prediction was determined by receiver operating characteristic (ROC) curves, survival curves, and nomograms. Finally, the correlation analysis was performed to investigate the relationships between the signature components of the ceRNA network and the immune cell signature.
RESULTS: Two completed survival analyses included selected components of the ceRNA network (ELN, SREBF1, DSC2, TTLL7, DIP2C, SATB1, hsa-miR-20a-5p, and hsa-miR-29c-3p) and selected immune cell types (M0 macrophages, M2 macrophages, resting mast cells, and neutrophils). ROC curves, survival curves (all P values <0.05), nomograms, and calibration curves indicated that the accuracy of the two survival analyses was acceptable. Moreover, the correlations between TTLL7 and resting mast cells (R=0.24, P<0.001), DSC2 and resting mast cells (R=-0.23, P<0.001), ELN and resting mast cells (R=0.44, P<0.001), and hsa-miR-29c-3p and M0 macrophages (R=-0.29, P<0.001) were significant, indicating that interactions of these factors may play significant roles in the prognosis of BC.
CONCLUSIONS: TTLL7, DSC2, ELN, hsa-miR-29c-3p, resting mast cells, and M0 macrophages may play an important role in the development of BC. However, additional studies are needed to confirm this hypothesis. 2021 Translational Andrology and Urology. All rights reserved.

Introduction

Bladder cancer (BC) is the ninth most common malignant tumor worldwide (1), accounting for an estimate of 549,000 new BC cases and 200,000 BC-related deaths in 2018 (2). BC is classified into non-muscle invasive bladder cancer (NMIBC) and muscle invasive bladder cancer (MIBC), and transitional cell carcinoma is the most common pathological type (3,4). NMIBC accounts for 75% of BC cases, and 50% of NMIBC cases progress to MIBC (5). Transurethral resection of bladder tumor (TURBt) is considered the main treatment for NMIBC, and radical cystectomy is the strategy for MIBC treatment (6). However, current treatment strategies for BC have not significantly improved the 5-year survival rate over the past decade (7). Moreover, BC has become one of the most expensive to treat solid tumors due to high recurrence rate (8). Therefore, it is essential to identify biomarkers that can accurately evaluate the diagnosis, treatment, and prognosis of BC.

Noncoding RNAs (ncRNAs), including long noncoding RNAs (lncRNAs) and microRNAs (miRNAs), have been recently recognized as the key regulatory factors in tumors (9). Overexpression of lncRNA DLX6 antisense RNA 1 (DLX6-AS1) and insulin-like growth factor binding protein 4-1 (IGFBP4-1) may lead to a poor prognosis in BC patients (10,11). Upregulation of miRNA-133b and miR-375-3p inhibits the proliferation and metastasis of BC cells (12,13). The occurrence of tumors is closely related to interactions between lncRNAs, mRNAs, and miRNAs (14). LncRNAs and mRNAs can competitively bind miRNAs by miRNA response elements (MREs) to generate a competitive endogenous RNA (ceRNA) network to regulate the expression of various RNAs and proteins (15,16). This ceRNA network can be used to investigate the functions of lncRNAs and to identify biomarkers related to the diagnosis and prognosis of BC, which is essential for assisting the clinicians with early diagnosis, risk assessment, and appropriate treatment decision-making for patients with BC (17,18).

Immunotherapy has become a promising antitumor strategy that recognizes tumors as foreign antigens and inhibits the proliferation and metastasis of the tumor cells by inducing active or passive immune responses (19,20). In BC, Bacillus Calmette-Guerin vaccine (BCG) suppresses tumor cells by activating immune cell infiltration and is the gold standard for therapy of high-risk NMIBC (21,22). However, approximately 40% of BC patients do not respond to BCG, and 15% of BC patients experience progression to MIBC after the treatment (23,24). Recent studies have shown that immune-related genes were significantly associated with the prognosis of cancer. High expression of ENDRA was associated with lower survival rate in patients with advanced BC (25). Overexpression of RAC3 promoted the migration of cancer cells and was associated with a decrease in the recurrence-free survival in ER3-positive breast cancer (26). The exploration of immune checkpoint inhibitors (ICIs) facilitated great progress in immunotherapy. In phase II clinical trials, the neoadjuvant use of ICIs resulted in pathological complete responses in patients with MIBC (27).

A number of recent studies demonstrated important roles of the ceRNA networks and immune cell infiltration in evaluation of prognosis in patients with malignant tumors (28,29). Some studies suggested that a combination of a ceRNA network and immune cell infiltration has potential value in cancer research (30,31). Therefore, we constructed a ceRNA network for BC using the data downloaded from The Cancer Genome Atlas (TCGA) database. Two survival analyses were performed, including analysis based on selected genes included in the ceRNA network and analysis based on the data on immune cell infiltration estimated by the Cell-type Identification By Estimating Relative Subsets Of RNA Transcripts (CIBERSORT) algorithm. These two survival analyses were used to identify the key factors and immune cells that can accurately predict patients with BC, and these factors were combined for subsequent analysis, which is very important for early diagnosis and treatment of BC. Finally, the relationships between these two survival analyses were analyzed (). The article is presented in accordance with the STROBE reporting checklist (available at http://dx.doi.org/10.21037/tau-20-1250).

Figure 1

The flow diagram of the present study. BC, bladder cancer; TCGA, The Cancer Genome Atlas; DEGs, differentially expressed genes; ceRNA, competitive endogenous RNA; CIBERSORT, cell-type identification by estimating relative subsets of RNA transcripts.

Methods

Identification of differentially expressed genes (DEGs)

The data for 411 BC tissues and 19 normal tissues were downloaded from the TCGA dataset (The Cancer Genome Atlas) (https://portal.gdc.cancer.gov/) and used for mRNA and lncRNA analysis. The data for 418 BC tissues and 19 normal tissues were downloaded from the TCGA dataset and used for miRNA analysis. The clinical data of 409 patients with BC were obtained from TCGA (including age, gender, grade, stage, and TNM classification) (). The DESeq2 method was used to identify DEGs, including lncRNAs, miRNAs, and mRNAs. DESeq2 can improve the stability and interpretability of estimation fold change (FC) and dispersion by shrinkage estimators (32,33). The data corresponding to duplicate and abnormal samples (not primary tumors and abnormal samples) were deleted, and the remaining data were normalized by the trimmed mean of M-values (TMM) method implemented in edgeR and further transformed by the voom method provided in limma. TMM normalization is a simple and effective method for estimation of relative RNA levels using RNA sequencing (RNA-seq) data. Moreover, voom is faster and more convenient than existing methods of RNA-seq analysis and converts RNA-seq data into a form that can be analyzed using tools similar to those used for microarray analysis. The cutoff criteria were defined as |log2 FC| >1.0 and false discovery rate (FDR) adjusted P<0.05. Ethical approval was not required because the data used in the present study were obtained from public databases. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Table 1

Characteristics of the included BC patients obtained from the TCGA database

Basic information	Total (n=409)
Age	69 (median)
Gender
Female	106
Male	303
Grade
High	385
Low	21
Unknow	3
Stage
I & II	132
III & IV	275
Unknow	2
T classification
T1 & T2	124
T3 & T4	253
TX	1
Unknow	31
N classification
N0	237
N1 &N2 & N3	131
NX	36
Unknow	5
M classification
M0	194
M1	11
MX	202
Unknow	2

BC, bladder cancer; TCGA, the The Cancer Genome Atlas.

Construction of a ceRNA network

StarBase was introduced to predict the lncRNA–miRNA and miRNA–mRNA interactions (34). StarBase (http://starbase.sysu.edu.cn/) is an open-source database that can display extensive and complex RNA-RNA (miRNA-lncRNA, miRNA-pseudogene, miRNA-circRNA, and miRNA-mRNA) and protein-RNA interaction networks by analyzing millions binding sites of Argonaute (Ago) and other RNA-binding protein (RBP) identified by crosslinking immunoprecipitation and high-throughput sequencing (CLIP-seq) (35). Interactions with P values <0.05 for the hypergeometric and correlation tests were considered significant. We used Cytoscape software version 3.8.0 to display the ceRNA network.

Survival analysis based on the components of the ceRNA network

All components of the ceRNA network were incorporated into a univariate Cox regression model to identify the genes related to survival (P<0.05). The survival analysis of the components of the ceRNA network was performed by least absolute shrinkage and selection operator (LASSO) Cox regression to avoid overfitting (36). The median value of the risk score calculated by LASSO was used as the cutoff value to divide patients into the high-risk and the low-risk groups. Kaplan-Meier survival curves were compared using the log-rank test and a nomogram to assess the survival outcomes [overall survival (OS)]. A receiver operating characteristic (ROC) curve and a calibration curve were constructed to evaluate the accuracy of the risk score.

Immune cell infiltration

We used CIBERSORT to evaluate the immune infiltration data for each BC sample downloaded from TCGA. CIBERSORT is a deconvolution algorithm that can predict the abundance of 22 immune cell types based on the gene expression data (37,38).

Survival analysis based on selected immune cells

All immune cells were incorporated into a univariate Cox regression model to identify immune cells related to survival (P<0.05). The survival analysis based on immune cells was performed by LASSO Cox regression to avoid overfitting. The patients were divided into the high-risk and low-risk groups based on the median value of the risk score. Kaplan-Meier survival curves were compared using the log-rank test and a nomogram to assess the survival outcomes (OS). A ROC curve and a calibration curve were constructed to evaluate the accuracy of the risk score. Finally, Pearson correlation analysis was used to assess the correlations between various immune cells and between the components of the immune cell signature and the components of the ceRNA network signature.

Multidimensional validation

Multidimensional verification was performed using a series of external databases, including Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/), the Gene Expression Profiling Interactive Analysis (GEPIA), and TargetScan database (http://www.targetscan.org/), to reduce possible errors. We also downloaded the relevant immunohistochemical images from the public database Human Protein Atlas database (HPA) (https://www.proteinatlas.org/) to detect the gene expression in tumor tissues and normal tissues, but no specific staining method was provided in the database.

Statistical analysis

All statistical analyses were performed using R 3.53 software. A ceRNA network was constructed using the GDCRNATools package of R (39). The log-rank test and univariate Cox regression analyses were performed using the survival package of R. The survival time and status of the patients were assumed to be dependent variables, and the expression of the components of the ceRNA network and immune cell infiltration were assumed to be independent variables in Cox regression used for two survival analyses. LASSO Cox regression was performed using the glmnet package of R. ROC curves were generated and plotted using the survivalROC package of R. Nomograms and calibration curves were generated using the rms package of R. The differences in the gene expression levels were evaluated by the Mann-Whitney U test.

Results

DEGs

We identified 2,581 differentially expressed mRNAs (1,261 upregulated and 1,320 downregulated mRNAs), 223 differentially expressed miRNAs (152 upregulated and 71 downregulated miRNAs), and 216 differentially expressed lncRNAs (39 upregulated and 177 downregulated lncRNAs) in BC versus normal tissues ().

Figure 2

DEGs. Differentially expressed mRNAs (A), lncRNAs (B), and miRNAs (C) were identified in BC. Red dots represent upregulated genes, and green dots represent downregulated genes. DEGs, differentially expressed genes; BC, bladder cancer.

Construction of a ceRNA network and survival analysis

We constructed a ceRNA network composed of 110 mRNAs, 15 miRNAs, and 7 lncRNAs (). Then, all components of the ceRNA network were incorporated into a univariate Cox regression model to identify components related to survival (P<0.05) (). The survival analysis by LASSO Cox regression used the selected components (; ). The median value of the risk score (1.072) calculated by LASSO was used as the cutoff value to divide BC patients into the high-risk and low-risk groups. ROC curve showed that OS of BC patients was perfectly predicted by the risk score [area under the curve (AUC) for 1-year survival =0.690; AUC for 3-year survival =0.707; and AUC for 5-year survival =0.743] (). The prognosis of BC patients in the high-risk group was worse than that in the low-risk group (P<0.001) (). A nomogram predicted the probability of 1-, 3-, and 5-year OS (). The calibration curve indicated that prediction using the nomogram was accurate ().

Figure 3

A ceRNA network. We constructed a ceRNA network using differentially expressed mRNAs, lncRNAs, and miRNAs. Purple ellipses represent mRNAs, green ellipses represent lncRNAs, and red ellipses represent miRNAs. The size of red ellipses corresponds to the number of the edges. CeRNA, competitive endogenous RNA.

Table 2

The genes related to survival identified by univariable Cox regression

Gene	HR	95% CI (low)	95% CI (high)	P
TIMP2	1.142	1.045	1.247	0.003
ELN	1.131	1.054	1.214	0.001
NTN1	1.099	1.024	1.179	0.009
MYLK	1.109	1.025	1.200	0.010
KIF26A	1.154	1.047	1.273	0.004
CYBRD1	1.120	1.015	1.236	0.024
SREBF1	1.174	1.021	1.349	0.024
TNS1	1.102	1.010	1.202	0.030
PTHLH	1.064	1.011	1.119	0.017
SEC23A	1.424	1.165	1.740	0.001
SRPX	1.129	1.047	1.218	0.002
RASD1	1.101	1.011	1.200	0.026
EFEMP1	1.119	1.056	1.186	0.000
AKT3	1.128	1.015	1.253	0.026
CTGF	1.121	1.032	1.217	0.007
RECK	1.147	1.021	1.288	0.021
PMEPA1	1.143	1.050	1.244	0.002
ATXN1	1.163	1.015	1.334	0.030
DSC2	1.113	1.041	1.191	0.002
PI15	1.100	1.023	1.183	0.010
THBS1	1.145	1.038	1.264	0.007
TTLL7	1.175	1.069	1.291	0.001
TPM1	1.185	1.060	1.325	0.003
ATP8B2	1.235	1.105	1.381	0.000
LRIG1	1.160	1.065	1.264	0.001
SETD7	1.252	1.019	1.538	0.032
DAAM2	1.139	1.022	1.269	0.019
LATS2	1.232	1.042	1.457	0.015
DIP2C	1.322	1.103	1.584	0.003
MSX1	1.162	1.050	1.285	0.004
FSTL1	1.192	1.053	1.350	0.005
EDIL3	1.121	1.013	1.241	0.027
ZCCHC24	1.135	1.018	1.265	0.023
JCAD	1.123	1.008	1.251	0.036
MAPRE2	1.181	1.025	1.360	0.021
ZEB2	1.124	1.003	1.259	0.044
PCDH7	1.068	1.001	1.139	0.045
SH3RF3	1.164	1.041	1.302	0.008
ZBTB4	1.349	1.037	1.754	0.026
SATB1	0.853	0.761	0.956	0.006
RFLNB	1.205	1.050	1.382	0.008
DUSP8	1.122	1.003	1.255	0.044
RUSC2	1.160	1.005	1.339	0.042
AC074117.1	0.746	0.607	0.918	0.006
hsa-let-7c-5p	1.161	1.067	1.264	0.001
hsa-miR-106b-5p	0.795	0.642	0.985	0.036
hsa-miR-17-5p	0.835	0.706	0.988	0.036
hsa-miR-20a-5p	0.804	0.688	0.939	0.006
hsa-miR-29c-3p	0.846	0.745	0.961	0.010
hsa-miR-93-5p	0.775	0.653	0.919	0.003

HR, hazard ratio; CI, confidence interval.

Figure 4

The ceRNA network signature. The ceRNA signature was constructed based on the components of the ceRNA network associated with survival according to LASSO Cox regression analysis (A,B). The results of ROC curve (C), Kaplan-Meier survival curve (D), nomogram (E), and calibration curve (F) analyses are shown. CeRNA, competitive endogenous RNA; LASSO, least absolute shrinkage and selection operator; ROC, receiver operating characteristic.

Table 3

The coefficients of included genes obtained from LASSO Cox regression

Gene	Coefficient
PARVB	0.0962
RAP1B	0.0529
PIK3CA	−0.0132
PGF	−0.0604
VEGFA	0.0294
SDC1	0.1073
SPP1	0.1006
FLNC	0.0940

LASSO, least absolute shrinkage and selection operator.

CIBERSORT was used to evaluate the immune infiltration data for each BC sample downloaded from TCGA, and bar plots and heatmaps were used to visualize the data (). The proportions of naive B cells (P=0.001), M0 macrophages (P<0.001), M1 macrophages (P<0.001), and resting mast cells (P=0.012) were different in the tumor and normal tissues (). Analysis of correlations between DEG expression and immune cell infiltration in BC and normal tissues () indicated that high levels of CD4 memory-activated T cells (P=0.033) and CD8 T cells (P=0.014) were associated with favorable OS, and high levels of memory B cells were associated with unfavorable survival (P=0.002) (). Additionally, we evaluated correlations between clinical and pathological factors and immune cells ().

Figure 5

Immune cell infiltration. The composition of immune cells was assessed by the CIBERSORT algorithm in BC (A,B). The proportions of naive B cells, M0 macrophages, M1 macrophages, and resting mast cells were different in BC and normal tissues (C). The results of the correlation analysis between DEG expression and immune cell infiltration are shown as a heatmap (D). High level of activated memory CD4 T cells (E) and CD8 T cells (F) was associated with favorable OS, and a high level of memory B cells was associated with unfavorable OS (G). CIBERSORT, cell-type identification by estimating relative subsets of RNA transcripts; BC, bladder cancer; OS, overall survival.

Figure 6

Clinical and pathological factors. The correlations between clinical and pathological factors and immune cells.

Survival analysis based on immune cells

All immune cells were incorporated into a univariate Cox regression model to identify immune cells associated with survival (P<0.05) (). LASSO Cox regression was used for survival analysis to avoid overfitting () (). The median value of the risk score (0.805) calculated by LASSO was used as the cutoff value to divide patients into the high-risk and low-risk groups. ROC curves indicated that OS of BC patients was perfectly predicted by the risk score (AUC for 1-year survival =0.704; AUC for 3-year survival =0.666; and AUC for 5-year survival =0.648) (). The prognosis of BC patients in the high-risk group was worse than that in the low-risk group (P=0.004) (). A nomogram predicted the survival probability of 1-, 3-, and 5-year OS (). The calibration curve confirmed that prediction using the nomogram was accurate (). Finally, correlations between the components of the ceRNA signature and immune cell signature were analyzed (). Tubulin tyrosine ligase like 7 (TTLL7) and resting mast cells (R=0.24, P<0.001), desmocollin 2 (DSC2) and resting mast cells (R=−0.23, P<0.001), elastin (ELN) and resting mast cells (R=0.44, P<0.001), and hsa-miR-29c-3p and M0 macrophages (R=−0.29, P<0.001) had significant correlations ().

Table 4

The immune cells related to survival identified by univariable Cox regression

Immune cell	HR	95% CI (low)	95% CI (high)	P
T cells CD8	0.041	0.004	0.477	0.011
T cells follicular helper	0.001	0.000	0.660	0.037
Macrophages M0	7.020	1.530	32.219	0.012
Macrophages M2	19.460	2.017	187.741	0.010
Mast cells resting	93.920	1.173	7,517.048	0.042
Neutrophils	38,449.382	6.225	237,470,827.229	0.018

HR, hazard ratio; CI, confidence interval.

Figure 7

The immune cell signature. The immune cell signature was generated by LASSO Cox regression analysis based on immune cells (A,B). The results of ROC curve (C), Kaplan-Meier survival curve (D), nomogram (E), and calibration curve (F) analyses are shown. LASSO, least absolute shrinkage and selection operator; ROC, receiver operating characteristic.

Table 5

The coefficients of included immune cells obtained from LASSO Cox regression

Immune cell	Coefficient	P
Macrophages M0	2.229	0.006
Macrophages M2	3.022	0.014
Mast cells resting	3.386	0.137
Neutrophils	13.294	0.004

LASSO, least absolute shrinkage and selection operator.

Figure 8

Coexpression analysis. The results of the correlation analysis of DEGs expression and immune cell infiltration are shown as a heatmap (A). The correlations between TTLL7 and resting mast cells (B), DSC2 and resting mast cells (C), ELN and resting mast cells (D), and hsa-miR-29c-3p and M0 macrophages (E) were significant. ELN, elastin; DSC2, desmocollin 2; TTLL7, tubulin tyrosine ligase like 7; DEGs, differentially expressed genes.

Multidimensional verification included a series of external databases (GEO, GEPIA, HPA, and TargetScan database) to verify the characteristics of significant factors in the coexpression test and survival analysis based on immune cells. Four GEO cohorts were used, including GSE13507, GSE7476, GSE32894, and GSE31684. The results obtained using the TCGA database indicated that the expression of ELN and TTLL7 was higher in normal tissues and the expression of DSC2 was higher in BC tissues. Analysis using the GEPIA database confirmed that both ELN and TTLL7 were expressed at a low level in tumor tissues, and the expression of DSC2 was upregulated. Analysis of the GSE13507 cohort indicated a lack of significant differences in the expression of ELN and DSC2, and the expression of TTLL7 was decreased in BC tissues compared with that in normal tissues. In the GSE7476 cohort, DSC2 was overexpressed and the expression of TTLL7 was decreased in BC tissues; however, the expression of ELN was not significantly different in tumor and normal tissues. The above results were summarized in . HPA database search results showed that the expression levels of ELN, DSC2, and TTLL7 were significantly increased in BC tissues (). Analysis of the data of multiple databases (TCGA and GSE31894) indicated that high levels of ELN, DSC2, and TTLL7 were significantly associated with poor prognosis of BC patients (). Finally, we used GSE31684 to validate the survival analysis based on immune cells. The prognosis of the high-risk group was worse than that of the low-risk group (P<0.001) (). A nomogram was able to predict the survival probability of 1-, 3-, and 5-year OS (). The calibration curve indicated that prediction using the nomogram was accurate (). Unfortunately, we were unable to verify the survival analysis based on the components of the ceRNA network due to a lack of databases that include mRNA, miRNA, lncRNA, and complete clinical data, which were present only in the TCGA database. Query of the TargetScan database identified binding sites for hsa-miR-29c-3p in ELN and DSC2, and a binding site for hsa-miR-374b-5p was identified in TTLL7 in the ceRNA network ().

Table 6

Multidimensional validation

Database	Gene	Tumor vs. normal	P
TCGA	ELN	Down regulated	<0.001
	DSC2	Up regulated	<0.001
	TTLL7	Down regulated	<0.001
GEPIA	ELN	Down regulated	<0.001
	DSC2	Up regulated	<0.001
	TTLL7	Down regulated	<0.001
GSE13507 (GEO)	ELN	NS	0.184
	DSC2	NS	0.687
	TTLL7	Down regulated	<0.001
GSE7476 (GEO)	ELN	Down regulated	0.009
	DSC2	Up regulated	0.036
	TTLL7	Down regulated	0.009

TCGA, the The Cancer Genome Atlas; GEPIA, Gene Expression Profiling Interactive Analysis; GEO, Gene Expression Omnibus; NS, no significance.

Figure 9

HPA. HPA database search results showed that the expression levels of ELN (A,B), DSC2 (C,D), and TTLL7 (E,F) were significantly increased in BC tissues. HPA, Human Protein Atlas database; ELN, elastin; DSC2, desmocollin 2; TTLL7, tubulin tyrosine ligase like 7; BC, bladder cancer.

Figure 10

Analysis of the data of multiple databases (TCGA and GSE31894). High levels of ELN, DSC2, and TTLL7 were significantly associated with poor prognosis of BC patients (A-F). Finally, we used GSE31684 to validate the survival analysis based on immune cells. The prognosis of the high-risk group was worse than that of the low-risk group (G). A nomogram was able to predict the survival probability of 1-year, 3-year, and 5-year OS (H). The calibration curve indicated that prediction using the nomogram was accurate (I). Query of the TargetScan database identified binding sites for hsa-miR-29c-3p in ELN and DSC2, and a binding site for hsa-miR-374b-5p was identified in TTLL7 in the ceRNA network (J). ELN, elastin; DSC2, desmocollin 2; TTLL7, tubulin tyrosine ligase like 7; BC, bladder cancer; OS, overall survival.

Discussion

The number of estimated new BC cases reported in 2019 was 80,470, and 17,670 BC-related deaths were registered. BC is characterized by complex biological behavior, a high rate of metastasis, and a high recurrence rate (40,41). After diagnosis, the physical, mental, and social health-related quality of life (HRQoL) of patients with BC significantly declines (42). Increasing evidence indicated that ceRNAs and dysregulation of the immune system are involved in the progression of human tumors (43-48).

Two survival analyses were performed in the present study based on the selected components of a ceRNA network (including ELN, SREBF1, DSC2, TTLL7, DIP2C, SATB1, hsa-miR-20a-5p, and hsa-miR-29c-3p) and selected immune cells (including M0 macrophages, M2 macrophages, resting mast cells, and neutrophils). Overexpression of SREBF1 was associated with a poor prognosis and promoted metastasis of esophageal carcinoma cells, and inhibition of SREBF1 augmented the efficacy of immune checkpoint blockades (49,50). Downregulated expression of DIP2C was detected in breast cancer, especially in basal-like and HER-2 subtypes (51). Numerous studies confirmed that SATB1 is involved in the occurrence and development of various human tumors (52-55). Hsa-miR-20a-5p inhibits epithelial mesenchymal transition (EMT) and invasion of tumor cells by targeting STAT3 in endometrial cancer (56). The infiltration of M2 macrophages expressing the CD163 antigen decreases the disease-free survival (DFS) rate in high-grade oral tongue squamous cell carcinoma and increases the mortality rate of patients with prostate cancer (57,58). The interaction between neutrophils and cancer is complex. Neutrophil-mediated tumor inflammation can promote the proliferation, invasion, and angiogenesis of tumor cells; however, some studies demonstrated that compounds produced by neutrophils can suppress tumor development (59-61).

The present study demonstrated significant correlations between TTLL7 and resting mast cells, DSC2 and resting mast cells, ELN and resting mast cells, and hsa-miR-29c-3p and M0 macrophages. Another study suggested that DS2 can be used as a new immunohistochemical biomarker for urothelial carcinoma (UC) with squamous differentiation; this cancer is more advanced and has worse prognosis than those in UC without squamous differentiation (62). In breast cancer, prostate cancer, and lung cancer, increased levels of mast cell infiltration were associated with improved survival rate of patients (63). Decreased expression of hsa-miR-29c-3p was associated with advanced clinical and pathological factors and poor prognosis in laryngeal squamous cell carcinoma, and infiltration of M0 macrophages was associated with poor OS, suggesting a negative correlation between these factors in tumors (64,65). Unfortunately, other components of ceRNA signature were rarely investigated in tumors, and previous studies did not assess possible relationships of these components with immune cells.

The present study has certain limitation. First, our data were obtained from public databases, and clinical and pathological factors included in these databases were not comprehensive; the number of tumor tissue samples was different from the number of normal tissue samples. Second, specific mechanism of infiltration of various immune cell types was not studied in detail. Third, the correlations between the components of the ceRNA network signature and the immune infiltrating cell signature were not confirmed in experiments.

Thus, two survival analyses performed in the present study were based on the selected components of the ceRNA network and selected immune cell types and were used to predict the prognosis of patients with BC; the nomograms containing both of these signatures may provide assistance to clinicians to improve individual management of BC patients. Moreover, significant correlations between TTLL7 and resting mast cells, DSC2 and resting mast cells, ELN and resting mast cells, and hsa-miR-29c-3p and M0 macrophages were detected, suggesting that these correlations play significant roles in the prognosis of BC. The present study provides some useful information for prediction of the prognosis of BC patients; however, further studies are needed to clarify the relationships between ceRNAs and immune infiltrating cells and relevant molecular mechanisms.

The article’s supplementary files as