Identification of Differentially Expressed irlncRNAs (DEirlncRNAs)
The flow chart of the study was shown in Fig. 1. Firstly, we downloaded transcriptome profiling consisting of 1053 breast cancer samples and 111 normal samples from the TCGA and then annotated according to gene transfer format (GTF) files from Ensembl. Next, we obtained the irlncRNAs by co-expression of immune-related genes from the ImmPort database. A total of 1041 irlncRNAs were identified (shown in Supplementary Table S1), and 55 were distinguished as DEirlncRNAs (Fig. 2a). Among the 55 distinguished DEirlncRNAs, 39 were upregulated while 16 were downregulated (Fig. 2b, Supplementary Table S2).
Construction Of Deirlncrna Pairs And A Prognostic Model
According to the iteration loop and the 0 or 1 matrix, we obtained 1172 valid DEirlncRNA pairs. Following, the DEirlncRNA pairs were merged with the survival information. Through Cox regression analysis, we obtained 20 DEirlncRNA pairs that were correlated with survival (Fig. 2c). Furthermore, 15 DEirlncRNA pairs were included in the model after Lasso regression analysis (Fig. 2d and e, Table 1).
Table 1
The risk model of 15 DEirlncRNAs pairs for breast cancer.
Pairs | Coef |
AP005131.7|LINC00511 | -0.01085 |
AP005131.7|AC009093.1 | -0.25014 |
U62317.4|AP000251.1 | -0.27626 |
AC011247.1|ATP2A1-AS1 | -0.17886 |
C6orf99|LINC01087 | 0.104848 |
C6orf99|U62317.1 | 0.294303 |
C6orf99|AC020663.2 | 0.272957 |
C6orf99|LINC02544 | 0.411963 |
LINC01929|AC020663.2 | 0.083829 |
LINC01929|AP005131.2 | 0.213263 |
U62317.1|ZNF350-AS1 | -0.05333 |
LINC00511|AP005233.2 | 0.010714 |
LINC00511|LINC01152 | 0.092366 |
AL356417.2|AP005131.2 | 0.337368 |
AC009093.1|AP005131.2 | 0.063 |
In order to verify the accuracy of the prognostic model, we drew the ROC curve and calculated the AUC value. We analyzed the model in the training set, and the 1-year AUC was 0.838, the 5-year AUC was 0.688, and the 10-year AUC was 0.788 (Fig. 3a). We validated the model in the validation set, and observed that all the AUC values were more than 0.680 (1-year AUC = 0.760, 5-year AUC = 0.719, 10-year AUC = 0.683, Fig. 3b). Similar trend was obtained for the total set (1-year AUC = 0.797, 5-year AUC = 0.700, 10-year AUC = 0.736, Fig. 3c).
According to the median value of risk score in the training set, all the samples were divided into the high-risk group and the low-risk group. We then performed a survival analysis of the training, the validation, and the total sets using Kaplan-Meier curves. As shown in Fig. 3d-f, patients in the high-risk group had worse overall survival than those in the low-risk group (P < 0.001).
Assessment of the Correlation between the Risk Model and Clinicopathological Characteristics
The risk curves and scatter plots were used to display the risk score and the survival outcome of each breast cancer patient in the training, the validation, and the total sets. The results showed that the mortality in the low-risk group was lower than in the high-risk group (Fig. 4a-f). The heatmap was used to show the relationship between DEirlncRNAs pairs and risk score in the three sets. AP005131.7 | LINC00511, AP005131.7 | AC009093.1, U62317.4 | AP000251.1, AC011247.1 | ATP2A1-AS1, U62317.1 | ZNF350-AS1 pairs had more values of 1 in the high-risk group than in the low-risk group. On the contrary, C6orf99 | LINC01087, C6orf99 | U62317.1, C6orf99 | AC020663.2, C6orf99 | LINC02544, LINC01929 | AC020663.2, LINC01929 | AP005131.2, LINC00511 | AP005233.2, LINC00511 | LINC01152, AL356417.2 | AP005131.2, AC009093.1 | AP005131.2 pairs had more values of 1 in the low-risk group than in the high-risk group (Fig. 4g-i).
To evaluate whether the risk model of the DEirlncRNAs pairs was an independent prognostic factor for breast cancer, univariate and multivariate Cox regression analyses were conducted. In the training set, the hazard ratio (HR) of risk score and 95% CI were 1.737 and 1.419–2.116 in univariate Cox regression analysis (P < 0.001), and 1.714 and 1.338–2.197 in multivariate Cox regression analysis (P < 0.001) respectively (Fig. 5a and b). This suggested that the risk model of DEirlncRNAs pairs were independent prognostic factors in patients with breast cancer. In the validation set, the risk model also showed statistical differences by univariate Cox regression analysis (P < 0.001, HR = 1.728, 95% CI [1.334–2.239]) and multivariate Cox regression analysis (P < 0.001, HR = 1.816, 95% CI [1.369–2.408], Fig. 5c and d). Similar results were obtained in the total set by univariate Cox regression analysis (P < 0.001, HR = 1.722, 95% CI [1.471–2.016]) and multivariate Cox regression analysis (P < 0.001, HR = 1.697, 95% CI [1.416–2.035], Fig. 5e and f).
Furthermore, we compared the differences in clinicopathological characteristics between the low-and high-risk group of all the samples. As shown by the strip chart and scatter plots, ER, PR, age, clinical stage, T stage, M stage, and N stage were significantly related to the risk score (Fig. 6).
Estimation of Tumor-Infiltrating Immune Cells and ICIs with Risk Assessment Model
Because lncRNAs play an important role in tumor immunity, we investigated the relationship between the model and the tumor immune microenvironment using all the breast cancer samples. As shown in Fig. 7a and b the immune score and stroma score were higher in the low-risk group than in the high-risk group (P < 0.05). We compared various tumor-infiltrating immune cells between the high-risk and low-risk groups, the results are listed in Supplementary Fig. S1. In addition, we integrated the relationships between various immune cells within the risk group performed by Spearman correlation analysis into a bubble chart. The results showed that the risk score was more positively associated with M2 macrophages and cancer-associated fibroblasts, whereas they were negatively associated with T cells, CD8 + T cells, NK cells, B cells, and M1 macrophages (Fig. 7c, Supplementary Table S3).
Immunotherapy is a novel and effective treatment method for breast cancer. We analyzed differences in ICIs common in breast cancer between high- and low-risk groups. The expression of PCDC1 (PD-1), CD274 (PD-L1), CTLA4, and CDK4 were significantly lower in the high-risk group (Fig. 8a-d).
Analysis of the Correlation between the Risk Model and Chemotherapeutic Drugs
The prediction of how effective a chemotherapy drug can guide the selection of clinical drugs. We assessed the correlation between the risk model and the efficacy of common chemotherapeutics used in the treatment of breast cancer. A higher IC50 of chemotherapeutic agents, such as methotrexate, doxorubicin, and gemcitabine were associated with the high-risk group (P < 0.05, Fig. 8e-g). While paclitaxel had no significance with the model (P > 0.05, Fig. 8H). The data indicated that the model might predict the treatment response to chemotherapy agents.
Gene Set Enrichment Analysis
The GSEA results indicated that the high-risk score group had markedly positive correlations with 6 enrichment pathways and negative correlations with 23 enrichment pathways (Supplementary Table S4). As shown in Fig. 9, the enriched KEGG pathways contained “glycolysis gluconeogenesis”, “citrate cycle TCA cycle”, “fructose and mannose metabolism”, “pentose phosphate pathway”, “terpenoid backbone biosynthesis” and “steroid biosynthesis”.