Using immune-related lncRNAs to construct novel biomarkers and investigate the immune landscape of breast cancer.

BACKGROUND: The role of immune-related long noncoding RNAs (irlncRNAs) in breast cancer (BRCA) is still unclear. Recently, studies have performed analyses based on the expression of irlncRNAs, however, in the present study, we used a novel method that did not require the specific expression levels of lncRNAs of BRCA patients.
METHODS: We downloaded transcriptome and clinical data of BRCA patients from The Cancer Genome Atlas (TCGA), obtained immune genes from the Immport database, and extracted immune genes and lncRNAs for correlation analysis. Then, the differential expression of irlncRNA pairs (IRLPs) was determined and the prognostic signature was established by the IRLPs. The immune cell abundance of the TCGA-BRCA cohort was downloaded from the Tumor IMmune Estimation Resource (TIMER) database, and the relationship between the risk score of the IRLP signature and immune cell abundance was analyzed. Finally, we explored the relationship between risk scores and drug sensitivity based on the R package pRRophetic.
RESULTS: Univariate cox regression results showed that 33 IRLPs had significant effects on the overall survival (OS) of BRCA patients. Then 22 IRLPs were obtained via lasso regression for further analysis. Multivariate regression analysis obtained 12 IRLPs to establish the IRLP prognostic signature. The model showed that this IRLP signature could act as a prognostic biomarker for BRCA patients. Kaplan-Meier (KM) survival analysis indicated that low-risk patients of IRLP's signature had a better OS (P<0.001). Advanced status BRCA patients may have higher risk scores, and univariate and multivariate cox regression analyses showed that risk scores were independent prognostic factors of clinical features (P<0.001). The results of the relationship between risk scores and immune infiltration showed that M1 macrophages were higher in the low-risk group (P=0.00015), while M2 macrophages were higher in the high-risk group (P=0.0015). The high-risk group had a greater sensitivity to chemotherapeutic agents such as cisplatin, docetaxel, doxorubicin, and gemcitabine.
CONCLUSIONS: In present study, we used a novel method that did not require the specific expression levels of lncRNAs of BRCA patients, which can be used as a novel model for predicting the prognosis of BRCA patients. 2021 Translational Cancer Research. All rights reserved.

Introduction

Breast cancer (BRCA) is the most common type of cancer in females worldwide. BRCA is highly heterogeneous, and there are significant differences in etiology, pathological manifestations, and the prognosis of individuals (1). Despite the developments of diagnostic methods and therapeutic regimens, the outcome of BRCA is still unsatisfactory, especially for progressed patients (2). Apart from being highly heterogeneous, the main reasons for this outcome may include metastasis, therapy resistance, and immune inhibition (3-5). Therefore, it is urgent to further understand the underlying mechanisms of these processes in BRCA, which may unveil potential therapeutic methods for advanced and refractory BRCA patients to lengthen their lifespan.

Long noncoding RNAs (lncRNAs) are a type of RNA which are characterized by more than 200 nucleotides in length, but do not code proteins. Recently, the dysregulation of lncRNAs has been observed in multiple cancer types, including colorectal cancer, lung cancer, and bladder cancer (6-8). Furthermore, lncRNAs act as important regulators of the progression of various cancers types (6,9,10). For BRCA, Liang et al.’s study indicated that lncRNA BCRT1 can promote BRCA progression via inhibiting the miR-1303/PTBP3 axis (11). Shi et al. demonstrated that lncRNA DILA1 can induce tamoxifen resistance via decreasing cyclin D1 degradation (12). However, Kim et al.’s study showed that lncRNA MALAT1 can inhibit BRCA metastasis (13). From these findings, we can conclude that lncRNAs may have dual roles in cancer processes.

To date, a number of studies have proven that the tumor microenvironment (TME) plays a significant role in BRCA progression (14-16). Apart from immune cells and stromal cells of the TME, several regulators such as lncRNAs have essential roles in the TME (17). Xu et al. showed that lncRNA SATB2-AS1 could regulate the tumor metastasis in colorectal cancer (18). LncRNA NEAT1 accelerates the progression of aggressive endometrial cancer by miR-361-regulated networks and TME-related genes (19). From these studies, lncRNAs may directly or indirectly affect the TME to regulate tumor processes. LncRNAs that regulate immune cell infiltration have been called immune-related lncRNAs (irlncRNAs). Several studies have proven that some irlncRNAs can be observed in BRCA patient samples such as lncRNA OSTN-AS1 and lncRNA TCL6 (20,21), and both of them affect the TME via regulating immune cell infiltration. Additionally, lncRNA TCL6 led to poor outcomes of BRCA patients (21), but the relationship between immune cell infiltration, lncRNAs, and drug resistance was not investigated. Compared with these studies, we explored the association between irlncRNAs and immune cell infiltration. We also constructed a novel modeling algorithm, estimated the predictive value, and explored drug sensitivity, as well as diagnostic sensitivity and specificity for BRCA patients.

We present the following article in accordance with the REMARK reporting checklist (available at https://dx.doi.org/10.21037/tcr-21-783).

Methods

Raw data

Transcriptome data and clinical data of 1,109 BRCA samples and 113 normal samples were extracted from The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/). An immune-related gene list was collected from the Immport website (https://www.immport.org/). Immune cell abundance of the TCGA-BRCA cohort was downloaded from the Tumor IMmune Estimation Resource database (TIMER, https://cistrome.shinyapps.io/timer/).

Obtaining irlncRNAs

We extracted the lncRNAs of the TCGA-BRCA cohort from the RNAseq expression matrix, and analyzed the correlation between immune genes and lncRNAs. The threshold values of correlation analysis were set as coefficient R >0.4 and P value <0.001 (Spearman).

Screening differentially expressed immune-related lncRNAs (irlncRNAs)

The differences in irlncRNAs were analyzed (tumor vs. normal) via the limma package in R software (https://www.r-project.org/). The filtering standards were log-fold change |logFC| >0.5 and false discovery rate (FDR) <0.05. The ggplot2 package was used to generate the heatmap plot and volcano plot.

Construction of immune-related lncRNA pairs (IRLPs)

The differentially expressed IRLPs were paired one by one. In each IRLP, the expression of the former lncRNA was higher than that of the latter lncRNA, then we recorded the expression of this IRLP as 1, or otherwise recorded as 0.

Construction of a prognostic signature for IRLPs

Firstly, we performed univariate cox regression on all IRLPs (screening P value <0.001), then we used lasso regression to obtain a more refined signature by constructing a penalty function. Finally, multivariate cox regression (stepwise method) was used to construct a prognostic signature based on the screening genes of lasso regression. The Kaplan-Meier (KM) method was used to generate the survival curve. The maximum Youden index (sensitivity + specificity-1) of the 1-year receiver operating characteristic (ROC) was identified as the best cutoff value based on the R packages survivalROC and survminer. A nomogram and calibration curve of the prognosis signature were constructed via the R package rms. We only included patients with complete survival data for signature construction. P<0.05 was considered statistically significant.

Correlation analysis between the risk score and immune cells

The correlation between the risk score of the prognostic signature and immune cells was analyzed based on the TCGA-BRCA cohort obtained from the TIMER database, and the correlation bubble chart between the risk score and immune cells was generated via ggplot2, ggtext, and scales packages (Spearman). P<0.05 was considered statistically significant.

Analysis of risk score and drug sensitivity

The R package pRRophetic is based on gene expression and drug sensitivity data from a large number of cancer cell lines, and these models are applied to gene expression data from primary tumor biopsies (22). We analyzed the relationship between the risk score of the signature and drug sensitivity based on the pRRophetic package. P<0.05 was considered statistically significant.

Statistical analysis

The differential irlncRNAs were screened through the R package limma (https://www.r-project.org/), and the thresholds were set as |logFC| >1 and FDR <0.05. The KM plot was used to assess the differential survival between the different groups, and log-rank P<0.05 was considered statistically significant. Univariate and multivariate (stepwise) cox regression were used to construct the prognostic signature. The Wilcoxon test and Spearman correlation were used for difference analysis and correlation analysis, respectively. P<0.05 was considered statistically significant.

Ethical statement

The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Results

Results of IRLPs

We extracted 13,413 lncRNAs from the TCGA-BRCA cohort to analyze the correlation with immune genes, than we obtained 420 irlncRNAs. The heatmap plot demonstrates the top 50 up-regulated and down-regulated genes (). Differential expression analysis indicated that there were 130 up-regulated lncRNAs and 45 down-regulated lncRNAs (). Subsequently, we paired 175 immune-related differentially expressed lncRNAs and obtained a total of 10,084 IRLPs.

Figure 1

The difference analysis of immune-related lncRNAs. Heatmap plot of the top 50 up-regulated and down-regulated lncRNAs, blue and red represent low and high expression, respectively (A); Volcano plot of differential lncRNAs, green and red represent down-regulated and up-regulated lncRNAs, respectively (B).

Prognostic signature results of IRLPs

Univariate cox regression results showed that 33 IRLPs had significant effects on the overall survival (OS) of BRCA patients. The results of lasso regression indicated that λ =−4.7 was the optimal value, and 22 IRLPs were obtained for further analysis (). Multivariate regression analysis obtained 12 lncRNA pairs to establish the IRLP prognostic signature ().

Figure 2

Prognostic signature of IRLPs. Lasso regression analysis indicated that λ =−4.7 (A) was the optimal value based on the 33 prognostic IRLPs (B); The forest map shows 12 IRLPs of the prognostic signature. Green and red respectively represent favorable IRLPs and harmful IRLPs to the overall survival of breast cancer patients (C). IRLPs, immune-related lncRNA pairs.

Optimal cutoff value determination and KM survival analysis of the IRLP signature

The 1-year ROC curve of the IRLP signature showed that the area under the curve (AUC) value was 0.832 and the optimal cutoff value of the risk score was 1.912 (). The AUC values of 1-, 3-, and 5-year ROC curves of the signature were 0.832, 0.803, and 0.794, respectively (). The risk score curve and survival status map showed that according to the optimal cutoff value, 1,090 BRCA patients were divided into a high-risk group (n=277) and low-risk group (n=813) (). KM survival analysis indicated that patients with low risk had a better OS for BRCA (log-rank P<0.001) ().

Figure 3

Optimal cutoff value determination and KM survival analysis of the IRLP signature. The 1-year ROC curve and the determination of the optimal cutoff value (A); the 1-, 3-, and 5-year ROC curves, the X-axis and Y-axis represent false positive rates and true positive rates, respectively (B); The risk score curve, the X-axis represents the increasing risk scores of patients, and the Y-axis represents risk scores (C); The survival status map, the X-axis represents the increasing risk scores of patients, the Y-axis represents survival time, and the green and red represent alive and dead, respectively (D); KM survival analysis of high- and low-risk patients of the IRLP signature (E). AUC, area under the curve; KM, Kaplan-Meier; IRLP, immune-related lncRNA pair; ROC, receiver operating characteristic.

Relationship between signature risk scores and clinical features

The summary heatmap of correlation analysis between risk scores and clinical features (). The results demonstrated that risk scores had no significant correlation with age (); and a higher risk scores were correlated with higher stage (P<0.001, ), tumor (T) size (P<0.001, ), and node (N) staging (P<0.05, ) of BRCA patients. Nevertheless, risk scores had no significant correlated with distant metastasis ().

Figure 4

Relationship between signature risk scores and clinical features. Summary heatmap of difference analysis of clinical features and risk scores (A); correlation analysis of risk scores with age (B), stage (C), tumor size (D), node staging (E), and distant metastasis (F) of breast cancer patients.

Prognostic analysis of signature risk scores and clinical features

Univariate and multivariate cox regression analyses of risk scores and clinical features indicated that risk scores were independent prognostic factors of clinical features (P<0.001, ). Compared with other clinical features of BRCA, the 1-year AUC value of the ROC curve of risk scores was 0.832, which was higher than T (0.729), N (0.654), M (0.580), stage (0.732), and age (0.794) features (). Risk scores were better prognostic factors for BRCA.

Figure 5

Prognostic analysis of signature risk scores and clinical features. Forest map of univariate cox regression analysis of risk scores and clinical features (A); Forest map of multivariate cox regression analysis of risk scores and clinical features (B); ROC curve of risk scores and clinical features, the X-axis and Y-axis represent false positive rates and true positive rates, respectively (C). ROC, receiver operating characteristic.

Relationship between risk scores and immune cells in BRCA

We downloaded the results of 7 types of software which contained various kinds of immune cells in BRCA from the TIMER database, and analyzed the correlation between risk scores and immune cells. The bubble chart demonstrated the correlation analysis results of significant differences between risk scores and immune cells ().

Figure 6

Relationship between risk scores and immune cells based on the TIMER database. Bubble chart of the X-axis represents the correlation R value, the right side of 0 represents a positive correlation, and the left side of 0 represents a negative correlation. The different colors of the Y-axis represent the calculation results of different software. TIMER, Tumor IMmune Estimation Resource (https://cistrome.shinyapps.io/timer/).

Based on the optimal cutoff value, risk scores were divided into high- and low-risk groups, and the relationship between high and low risk and M1 macrophages and M2 macrophages was analyzed. The results showed that the M1 macrophage content in BRCA was higher in the low-risk group (P=0.00015, ), and M2 macrophage content in BRCA was higher in the high-risk group (P=0.0015, ).

Figure 7

Comparison of the content of macrophages between high- and low-risk groups. M1 macrophages (A); M2 macrophages (B).

Results of risk scores and drug sensitivity analysis

We analyzed the half maximal inhibitory concentration (IC50) of some drugs which are commonly used as chemotherapeutic agents in BRCA based on the R package pRRophetic. The results indicated that doxorubicin (P=0.0011, ), cisplatin (P=0.032, ), docetaxel (P=0.061, ), and gemcitabine (P=2e-07, ) had lower IC50 values in BRCA patients with low risk, suggesting that these drugs have better efficacy in the low-risk group.

Figure 8

Relationship between risk scores and drug sensitivity. Doxorubicin (A); cisplatin (B); docetaxel (C); and gemcitabine (D).

Discussion

With the better understanding of the association between lncRNAs and the immune system (23,24), more and more evidence indicates that lncRNAs can not only act as diagnostic biomarkers, but also prognostic signatures for cancer patients. Several studies have focused on using coding genes and noncoding RNAs to construct models to assess the prognostic and diagnostic values for BRCA (25-27). Most of the signatures in these studies were determined via quantifying transcriptional expression levels. In the present study, instead of using transcriptional expression levels, we used immune-related gene pairing to generate a reliable model with lncRNA pairs.

We downloaded transcriptome and clinical data of BRCA patients from TCGA, obtained immune genes from the Immport database, and extracted immune genes and lncRNAs for correlation analysis. Then, the differential expression of IRLPs was determined and the prognostic signature was established by the IRLPs. The immune cell abundance of the TCGA-BRCA cohort was downloaded from the TIMER database to analyze the relationship between the risk score of the IRLP signature and immune cell abundance. Finally, we explored the relationship between risk scores and drug sensitivity based on the R package pRRophetic. Shen et al. used the transcriptional expression levels of 11 irlncRNAs as a biomarker for predicting the prognosis of BRCA patients (28). It is now known that lncRNAs are the most abundant non-coding RNAs, with a number of significant biological functions. Our study indicated that we could select significant IRLPs from the initial identification differentially expressed irlncRNAs (DEirlncRNA). We could detect significant differential expression of IRLPs but not every lncRNA. Clinicians can take advantage of this novel model to divide patients into low or high risk for BRCA. Since these lncRNAs are related to immune-related genes, they may have a deep influence on the TME of BRCA. In our study, we found that the IRLPs TMPO−AS1|SIAH2−AS1, AP000695.1|AP005131.2, AL121832.2|AC109322.1, AC007637.1|ROCR, LINC01235|U62317.4, LINC01235|ROCR, LINC01235|MAGI2−AS3, LINC01235|U62317.1, AC005041.3|AP005131.2, CYTOR|MIR200CHG, AC004585.1|AL645608.7, and AP005131.2|AC245014.3 had significant roles in the prognosis of BRCA patients via univariate analysis. We further performed multivariate analysis which showed that TMPO−AS1|SIAH2−AS1, AL121832.2|AC109322.1, AC007637.1|ROCR, LINC01235|MAGI2−AS3, CYTOR|MIR200CHG, and AC004585.1|AL645608.7 had essential roles in the prognosis of BRCA patients. Several studies have proven that lncRNA TMPO−AS1 is a cancer promoter for BRCA (29,30). Additionally, lncRNA TMPO−AS1 plays a key role in ovarian cancer and bladder cancer (31,32). MAGI2−AS3 is a tumor suppressor for BRCA, and Du et al. showed that MAGI2−AS3 inhibited the migration and invasion of BRCA cells by sponging microRNA-374a (33). Liu et al. showed that CYTOR can enhance the tamoxifen resistance of BRCA cells through sponging miR-125a-5p (34). Besides having a significant effect on BRCA patients, Moradi et al. showed that CYTOR also acted as a diagnostic biomarker (35). Furthermore, a review by Liang et al. showed that CYTOR was an adverse signature for multiple cancer types (36). Other irlncRNAs have not been investigated in regards to their biological functions in cancer.

We further used lasso regressionto increase the predictive model accuracy and efficacy. We also improved the procedure of modeling by calculating every AUC value to identify the maximum value for an optimal model followed by comparison with other clinical parameters. The AUC values were also used to achieve an optimal point for model fitting instead of distinguishing the risk just by the median value. Apart from exploring the relationship between risk scores and the prognosis of BRCA patients, we also reanalyzed the outcome of BRCA patients based on clinical features such as stage, age, and gender via univariate and multivariate analyses. We also analyzed the relationship between immune cell infiltration and risk scores via multiple software. The final goal of cancer research is always to find novel potential therapeutic regimens for cancer patients, so we further investigated the drug sensitivity of BRCA. The results showed that DEirlncRNA pairs were more positively related to tumor-infiltrating immune cells, and have a positive correlation with common lymphoid progenitors, neutrophils, CD4+ Th2 T cells, macrophages, activated NK cells, macrophage M0 infiltration, but a negative correlation with CD4+ effector memory T cells, CD4+ central memory T cells, CD4+ naive T cells, B cells, and activated myeloid dendritic cells, among others. From the various immune cells of the TME for BRCA, our past study showed that macrophages play a key role in BRCA progression (37). In the present study, we mainly focused on the role of macrophage cell types in BRCA. The results showed that the high-risk group had higher M2 macrophage infiltration than the low-risk group, while the low-risk group had higher M1 macrophage infiltration. These results indicated that M2 macrophages may act as suppressors of BRCA while M1 macrophages may act as promoters of BRCA. Chen et al. showed that M2 macrophages promoted BRCA metastasis by the M2 macrophage-produced CHI3L1 protein (38). Sousa et al. revealed that BRCA cells can induce macrophages toward M2 polarization (39). These results indicate that BRCA can adjust the macrophage status to affect the tumor process. Li et al. showed that promoting M1 macrophages could inhibit the growth of BRCA cells via activating the NF-κB pathway (40). However, Guo et al. showed that M1 macrophages can induce BRCA stem cells by regulation of the Lin-28B-let-7-HMGA2 axis (41). These findings indicate that the roles of macrophages may be adjusted according to different TMEs. Wang et al. indicated chemotherapy and immunotherapy therapy responses of cancer via immune scores. Furthermore, a number of studies have shown that the TME has a significant influence on therapy response (42-44). In our study, our model showed that the high-risk group was related to sensitivity to chemotherapeutic regimens such as cisplatin, docetaxel, doxorubicin, and gemcitabine. Our results recommend clinicians to select chemotherapeutic regimens based on risk scores to obtain better therapeutic responses.

Conclusions

In conclusion, the present study showed that a novel signature constructed by irlncRNAs, that did not require an analysis of lncRNA expression levels, could assess the prognosis of BRCA patients and help clinicians to select chemotherapeutic regimens based on risk scores.