The level 3 mRNA expression profiles of BC were obtained from The Cancer Genome Atlas (TCGA) by the “TCGAbiolinks” R package (16). A total of 1,072 BC samples were found with available lncRNA and mRNA expression profiles, and 1,066 BC samples were found with available miRNA expression profiles. The lncRNA and mRNA expression data were converted from “fragments per kilobase million” into “transcripts per kilobase million”. We gathered clinical data on 1,087 BC patients, including their overall survival and progression-free survival data.

The single-sample Gene Set Enrichment Analysis (ssGSEA) model was used to calculate the amounts of immune cells and to convert the genomic data into the values of 28 different categories of human immune cells (17), which included the cells that were thought to be anti-tumor such as “activated CD4 T cell” and “activated CD8 T cell”. Several pro-tumor immunity cells, including “regulatory T cell” and “type-2 T-helper cell”, were also obtained. Two groups of BC patients were established according to the median value of “activated CD8 T cell”. The microenvironment cell populations counter (MCP-counter) method (18), which allows the accurate quantification of immune and stromal cells by genomic data, was used in this study. The ESTIMATE, a tool for predicting tumor purity and the presence of infiltrating stromal/immune cells in tumor tissues using gene expression data, was used (19).

The expression profiles of 8,633 lncRNAs were obtained using the annotation profile from the GENCODE program. After excluding 25% of lncRNAs having the lowest mean expression value, 6,475 lncRNAs were retained for further study. The expression profiles of 2,236 miRNAs were obtained by the miRNA annotation from the miRBaseVersions.db R package. After removing 25% of lncRNAs with the lowest mean expression value, 1,677 miRNAs were retained for further investigation. The expression profiles of 19,618 mRNAs were obtained by the mRNA annotation of the GRCh38 project. After excluding 25% of mRNAs with the lowest mean expression value, 14,714 mRNAs were retained for further study.

The differentially expressed RNAs (DElncRNAs, DEmiRNAs, and DEmRNAs) between two T-cell groups were calculated by the R package edgeR (20). DElncRNAs and DEmiRNAs were filtered using p-value of 0.05 and |log2fold change| > 0.3 as criteria. DEmRNAs were filtered using p-value of 0.05 and |log2fold change| > 0.5 as cutoffs. Upregulation mRNAs (UmRNAs), down-regulation mRNAs (DmRNAs), upregulation miRNAs (UmiRNAs), downregulation miRNAs (DmiRNAs), upregulation lncRNAs (UlncRNAs), and downregulation lncRNAs (DlncRNAs) were defined by their expression values in the high CD8 T-cell group.

Based on ceRNA theory, we constructed lncRNA–miRNA–mRNA ceRNA networks by lncRNA–miRNA and miRNA–mRNA correlations. First, we downloaded all the lncRNA–miRNA and miRNA–mRNA interactions that were downloaded from StarBase (http://starbase.sysu.edu.cn/) (21). The first ceRNA network that has a positive correlation value with CD8 T cells was constructed by selecting ncRNA–miRNA and miRNA–mRNA interactions that contained UmRNAs, DmiRNAs, and UlncRNAs. The second ceRNA network that has a negative correlation value with CD8 T cells was constructed by selecting ncRNA–miRNA and miRNA–mRNA interactions that contained DmRNAs, UmiRNAs, and DlncRNAs. In this study, Cytoscape software was utilized to visualize the network.

Since functional enrichment analysis is crucial for interpreting high-throughput omics data in life science, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed in R using the function of clusterProfiler (p-value <0.05) (22). The lncRNAs were ranked by the number of ncRNA–miRNA interactions. The top 4 lncRNAs were selected from the first ceRNA network, and the top 5 lncRNAs were selected from the second ceRNA network. These nine lncRNAs were defined as the hub lncRNAs.

Using the expression profiles of nine hub lncRNAs, machine learning models were built for the prediction of CD8 T cells. BC patients were equally and randomly divided into training and testing samples. Different machine learning models, including decision tree (DT), gradient boosting machines (GBM), generalized linear models (GLM), artificial neural networks (ANN), random forest (RF), and support vector machine (SVM), were chosen to build models to predict the CD8 T-cell group by the R “Caret” package (23). The parameters for each machine learning model were determined by the fivefold cross-validation. In order to evaluate the robustness and prediction accuracy of constructed models in BC, the “ROCR” package was used to calculate the value of the area under the ROC curve (AUC).

These nine hub lncRNAs were used to construct the Cox model. The coefficients for nine lncRNAs were obtained. The risk score of the samples was calculated by the coefficients and lncRNA expression profiles. Based on the median value, BC individuals were categorized into high- and low-risk classes. The correlation analysis of clinical information includes overall survival, age, gender, stage, and T, N, and M classifications in the TNM system. The log-rank test was used for the Kaplan–Meier curves of patient survival analyses.

Total RNA was extracted from patients’ breast cancer at the First Affiliated Hospital of Nanjing Medical University using TRIzol reagent (Invitrogen, CA, USA). Isolated RNA was reverse transcribed into cDNA using HiScript II Q RT SuperMix for qPCR (Vazyme Biotech Co.,Ltd. Nanjing, China) following standard protocols. Real-time quantitative PCR (qPCR) was performed with synthetic primers and ChamQ SYBR qPCR Master Mix (Vazyme Biotech Co., Ltd. Nanjing, China) with a Quant Studio 5 Real-Time PCR Detection System (Thermo Fisher Scientific, MA, USA). The relative expression levels of XIST were calculated and quantified with the 2^−ΔΔCt method after normalization with the reference β-actin expression. All primers used are listed in Table 1 .

The flowchart of this study is presented in Figure 1 . By the ssGSEA and MCP-counter methods, values of CD8 T cells were calculated. High values of CD8 T cells were correlated with better BC prognosis ( Supplementary Figure 1 ). Then, the BC samples were divided into two groups (low and high) by the median value of CD8 T cells from the ssGSEA algorithm. Based on the criteria of p-value <0.05 and |log₂FC| > 0.3, 1,599 DElncRNAs (811 downregulation and 788 upregulation, Figure 2A ) and 89 DEmiRNAs (11 downregulation and 78 upregulation, Figure 2B ) between low and high CD8 T-cell groups were obtained. Based on the criteria of p-value <0.05 and |log₂FC| > 0.5, 1,794 DEmRNAs (659 downregulation and 1,135 upregulation, Figure 2C ). Heatmaps of the top DElncRNAs, DEmiRNAs, and DEmRNAs are presented in Figures 2D–F . Upregulation mRNAs (UmRNAs), downregulation mRNAs (DmRNAs), upregulation miRNAs (UmiRNAs), downregulation miRNAs (DmiRNAs), upregulation lncRNAs (UlncRNAs), and downregulation lncRNAs (DlncRNAs) were defined by their expression values in the high CD8 T-cell group.

We downloaded the lncRNAs–miRNAs and miRNAs–mRNAs interactions from the StarBase database. Due to the ceRNA theory, lncRNAs and mRNAs have a positive regulatory interaction, while miRNAs and mRNAs have a negative regulatory interaction. Then, lncRNA–miRNA and miRNA–mRNA interactions that contain UmRNAs, DmiRNAs, and UlncRNAs were used to construct a CD8 T-cell positive ceRNA network ( Figure 3A ). The lncRNA–miRNA and miRNA–mRNA interactions that contain DmRNAs, UmiRNAs, and DlncRNAs were used to construct a CD8 T-cell negative ceRNA network ( Figure 3B ).

To outline the potential function of the genes in the ceRNA networks, the pathways from the CD8 T-cell positive ceRNA network and the CD8 T-cell negative ceRNA network were obtained. Top GO terms and KEGG pathways in the CD8 T-cell positive ceRNA network included many immune pathways such as “antigen processing and presentation” ( Supplementary Table 1 ). Top GO terms and KEGG pathways in the CD8 T-cell negative ceRNA network contained many cancer pathways such as “proteoglycans in cancer” ( Supplementary Table 2 ).

To predict the CD8 T-cell group of BC samples, we constructed a model based on the hub lncRNAs from the CD8 T-cell-positive/negative ceRNA networks. The top 4 lncRNAs from the CD8 T-cell-positive ceRNA network, and the top 5 lncRNAs from the CD8 T-cell-negative ceRNA network were selected. These nine lncRNAs (NEAT1, XIST, NORAD, MALAT1, MIR29B2CHG, LINC00943, AC005332.4, AC092718.4, and LINC01146) were defined as the hub lncRNAs, and the expression profiles of these hub lncRNAs were kept for further analysis. The expression profiles were evenly and randomly split into a training (50%) and a testing set (50%). Fivefold cross-validation was conducted on the training set to select the best parameters for DT, GBM, GLM, ANN, RF, and SVM. The models were trained on the training set with the selected best parameters. The log2CP value from DT was selected as “−9.65”, since it achieved the highest AUC value of 0.794 in the fivefold cross-validation ( Figure 4A ). The interaction depth value from GBM was selected as “7”, since it achieved the highest AUC value of 0.857 in the fivefold cross-validation ( Figure 4B ). The alpha value from GLM was selected as “0.7”, since it achieved the highest AUC value of 0.838 in the fivefold cross-validation ( Figure 4C ). The size value from ANN was selected as “11”, since it achieved the highest AUC value of 0.837 in the fivefold cross-validation ( Figure 4D ). The mtry value from RF was selected as “3”, since it achieved the highest AUC value of 0.842 in the fivefold cross-validation ( Figure 4E ). The log2C value from SVM was selected as “−1”, since it achieved the highest AUC value of 0.843 in the fivefold cross-validation ( Figure 4F ). Prediction accuracy was then measured by the AUC value calculated on the test set. On the testing set, the AUC values for DT, GBM, GLM, ANN, RF, and SVM were 0.742, 0.854, 0.848, 0.855, 0.852, and 0.843 ( Figures 5A–F ). The calculation process of the final model for DT is presented in Supplementary Figure 2 .

Thereafter, we developed a nine-lncRNA prognostic model for overall survival prediction of BC patients. The risk score = (0.61 * MIR29B2CHG) + (0.027 * NEAT1) + (0.006 * MALAT1) + (−1.21 * LINC00943) + (0.24 * LINC01146) + (2.1 * AC092718.4) + (−1.02 * AC005332.4) + (0.99 * NORAD) + (0.193 * XIST). The distribution of prognosis, including survival time and status in risk score groups, is presented ( Figure 6A ). Additionally, the increased expression of AC005332.4 and LINC00943 was observed in the low-risk BC samples. On the other hand, MIR29B2CHG, MALAT1, XIST, NORAD, and AC092718.4 expression levels were greater in the high-risk BC samples ( Figure 6B ). Low-risk individuals in the BC cohort survived longer than high-risk ones when it came to OS (p-value = 0.001, Figure 6C ). MIR29B2CHG, NEAT1, MALAT1, AC005332.4, NORAD, and XIST were elevated in normal breast samples than in breast cancer samples ( Figure 6D ). LINC01146 and AC092718.4 were more elevated in breast cancer samples than in normal breast samples ( Figure 6D ).

We analyzed the associations among nine lncRNAs (NEAT1, XIST, NORAD, MALAT1, MIR29B2CHG, LINC00943, AC005332.4, AC092718.4, and LINC01146), risk score, tumor purity, and immune cell infiltration ( Figure 7A ). NEAT1, XIST, NORAD, MALAT1, and MIR29B2CHG were significantly negatively related to CD8 T cells, cytotoxic lymphocytes, and T cells. These five lncRNAs were significantly positively associated with fibroblasts, endothelial cells, and neutrophils. On the other hand, LINC00943, AC005332.4, AC092718.4, and LINC01146 were significantly positively associated with most types of immune cells, including T and B cells. Risk score value was found to be negatively related to immune cells but positively associated with tumor purity. NEAT1, XIST, NORAD, MALAT1, MIR29B2CHG, and risk value were significantly negatively related to immune checkpoint genes including LAG3, TIGIT, CTLA4, and PDCD1 ( Figure 7B ). LINC00943, AC005332.4, AC092718.4, and LINC01146 were significantly positively associated with LAG3, TIGIT, CTLA4, and PDCD1.

Among nine lncRNAs, MIR29B2CHG, NEAT1, MALAT1, NORAD, and XIST were elevated in the low CD8 T-cell group than in the high CD8 T-cell group ( Supplementary Figure 3 ). LINC00943, LINC01146, AC005332.4, and AC092718.4 were elevated in high CD8 T cells group than in low CD8 T cells group ( Supplementary Figure 3 ). The associations between the risk model and clinical features were explored. Clinical characteristics, such as age, stage, T stage, N stage, and M stage, had no effect on the risk score ( Supplementary Figures 4A–E ). We found that a higher risk score was correlated with a negative prognosis on progression-free survival ( Supplementary Figure 4F ). Based on the best cutoff values of these nine lncRNAs, the samples were divided into high and low groups. The survival results showed that AC005332.4, AC092718.4, NORAD, and LINC00943 have significant correlations with prognosis ( Supplementary Figure 5 ).

To validate XIST expression in BC clinical samples, 10 paired BC and normal breast tissues were collected and compared using real-time quantitative PCR (RT-qPCR). As shown in Figure 7C , expression of XIST is significantly lower in BC samples compared to the adjacent normal breast tissues (p-value <0.05).