The transcriptome RNA-seq data [fragments per kilobase of transcript per Million mapped reads (FPKM)] of 1,222 samples were obtained from the TCGA public database, including 1110 BC tissues and 112 normal adjacent tissues (https://portal.gdc.cancer.gov/repository). The corresponding clinicopathological data of BC were also downloaded from TCGA together. The protein-coding genes and lncRNAs were distinguished by applying the ensembl human genome browser GRCh38.p13. A total of 170 NETs-related genes were acquired from previously published studies (as shown in Supplementary Table S1) (Dwyer et al., 2014; Papayannopoulos, 2018). Then the correlation between the expression of NETs-related genes and corresponding lncRNAs was quantified by calculating the Pearson correlation coefficients. The strict criteria were used to identify NETs-related lncRNAs, with p < 0.001 and the absolute value of Pearson correlation coefficient more than 0.4 (|R| > 0.4).

The expression profiles of NETs-related lncRNAs of the 112 normal breast samples and 1110 BC samples were obtained, and the differential expression analysis was carried out with |Log₂ fold change [FC]| >1 and false discovery rate (FDR) <0.05 using the “limma” and “pheatmap” R package. After deleting patients with incomplete information, 937 patients with breast cancer in total were randomly divided in an 8:2 ratio into training and validation cohorts for constructing and validating the NETs-related lncRNAs signature. Univariate Cox regression analysis for overall survival (OS) was employed to determine prognostic lncRNAs with p < 0.05 in the training cohort. The least absolute shrinkage and selection operator (LASSO) Cox regression algorithm was performed to lessen the chance of overfitting. Furthermore, multivariate Cox regression analysis was applied to calculate the regression coefficient of the prognostic risk score model.

The lncRNA-mRNA co-expression network was constructed to demonstrate the correlation between the NETs-related lncRNAs and their corresponding mRNAs, and visualized using the Cytoscape software (version 3.7.2, http://www.cytoscape.org/).

The risk score for each BC patient in both the training and validation cohorts was calculated according to the following equation: R i s k   s c o r e = ∑ i = 1 ∞ c o e f i × x i (1) where coefi is the regression coefficient and xi is the corresponding lncRNA expression level. The BC patients in the training cohort were separated into high- and low-risk groups based on the median risk score. The same median score was used to divide patients in the validation cohort into high- and low-risk groups. The Kaplan-Meier method, risk score heatmap, distribution of risk score and survival status were used to assess the validity of the prognostic risk model by applying the “survivalROC” and “pheatmap” R package.

Univariate and multivariate Cox regression algorithm were conducted to assess the prognostic significance of risk scores based on other clinical parameters (age, stage, and subtype) in the training cohort. The receiver operating characteristic (ROC) curves and C-index were performed to evaluate the predictive power of this signature. Then, a nomogram was constructed based on these parameters as implemented in the “rms” R package. The validity of prognostic model in their accuracy of prediction for 1-, 3-, and 5-years OS was further assessed by the calibration curves and time-dependent receiver operating characteristic (ROC) curves in the training and validation cohorts.

To investigate the potential biological functions of the 10 NETs-related lncRNAs in BC, the gene set enrichment analysis (GSEA) was conducted to look for the tumor hallmarks associated with risk scores. The FDR <0.25 and P adjusted value <0.05 were adopted as the criteria for statistical significance. An enrichment lot was performed to visualize the top five functions enriched by two groups. After that, the differentially expressed genes between the two risk scores groups were determined by differential expression analysis of two groups with |log FC| >1 and p < 0.05. The pathway enrichment analysis for Gene ontology (GO) and kyoto encyclopedia of genes and genomes (KEGG) on the screened genes were conducted by employing the “clusterProfiler” R package. Then, the “enrichplot” and “ggplot2” R packages were applied to visualize the enrichment results.

The 977 single nucleotide variants (SNVs) data were downloaded in the TCGA public database and the “Maftools” R package was applied to calculate the tumor mutation burden (TMB) for each BC patient. The correlation between the risk score and TMB was analyzed using the Spearman’s algorithm. The Kaplan-Meier survival analysis was performed for different TMB groups. Furthermore, the CIBERSORT deconvolution algorithm with 1,000 permutations was applied to calculate the abundance of immune cell infiltration in the tumor immune microenvironment of high- and low-risk patients for BC with p < 0.05.

Single-sample Gene Set Enrichment Analysis (ssGSEA) was performed to determine the biological function differences between the high- and low-risk groups by applying the “GSVA” R package. The potential immune checkpoint genes were chosen according to previous published literature. The different expression levels of immune checkpoint genes between risk groups were analyzed by the Wilcoxon test.

To assess the therapeutic efficacy of chemotherapy and targeted agents in constructing prognostic risk signature of NETs-related lncRNAs, the half inhibitory concentration (IC₅₀) of chemotherapeutic and target therapeutic drugs were calculated using the “ggplot2” and “pRRophetic” R packages. The Wilcoxon signed-rank method was applied to compare the IC₅₀ of different risk groups. Then, immunotherapy responsiveness was predicted in different risk groups by applying the Tumor Immune Dysfunction and Exclusion (TIDE) online Tool (http://tide.dfci.harvard.edu/).

A total of 10 matched pairs of tumor specimens and adjacent normal tissues were obtained from BC patients who underwent tumor resection. All tissue specimens were collected from the Breast Surgery Department of Sun Yat-sen University Cancer Center, Guangzhou, China, with patients’ consent and approval from the Medical Ethics Committee of the hospital.

For RNA expression assay, total RNA was extracted using TRIzol reagent (Invitrogen, 15596018) and quantified using a Nanodrop (Thermo Scientific). Then, cDNA was obtained using the PrimeScript™ RT Reagent Kit (Takara, RR036A) according to the manufacturer’s instruction. Quantitative real-time polymerase chain reaction (qRT-PCR) was performed with the TB Green™ Premix Ex Taq™ (TaKaRa, RR420A) using SYBR Green (Roche) on a LightCycler 480 (Roche). The relative abundance of each lnRNA, using GAPDH as an endogenous control, was calculated by 2^−∆∆CT method. The primers sequences used in this study are listed in Supplementary Table S2. Analysis between the two groups was performed by paired-sample t-tests. p < 0.05 was considered statistically significant.

A flow chart of the study is depicted in Figure 1. A total of 14,142 lncRNA transcripts and 19,658 protein-coding genes were identified from the TCGA database. The gene expression of NETs in BC was screened by matching the mRNA expression matrix of TCGA and 170 NETs-related genes. Afterwards, 655 lncRNAs were identified to be highly correlated with the expression of NETs-related genes (|R| > 0.4 and p < 0.001), among which 186 lncRNAs (as shown in Supplementary Table S3) were identified based on differential expression between normal and tumor tissue. A heatmap was performed to demonstrate the differential expression of the top 20 up-regulated and 20 down-regulated lncRNAs (Figure 2A), and volcano plot of differentially expressed lncRNAs showed 105 up-regulated lncRNAs and 81 down-regulated lncRNAs in BC (Figure 2B).

Next, to establish a convincing risk predictive model, BC patients were randomly divided in an 8:2 ratio into training (n = 749) and validation cohorts (n = 188). In the training cohort, a univariate Cox regression algorithm was performed to determine the lncRNAs related to BC patient survival and 18 prognosis-related candidate lncRNAs were filtered out (Figures 3A, B). The 18 lncRNAs were incorporated into the LASSO Cox regression algorithm to improve model accuracy and reduce model overfitting (Figure 3C), and cross-validation was performed and 13 prognostic lncRNAs were screened out (Figure 3D). To further improve clinical utility, 13 lncRNAs signatures were purified by performing stepwise multivariate Cox regression analysis (Figure 3E). Then, a prognostic risk model was constructed based on the correlation coefficient of the expression of lncRNAs in the multivariable Cox regression model. The following formula was applied to determine the risk score: risk score = (0.366 × SFTA1P expression level) + (−0.500 × ACTA2-AS1 expression level) + (−0.907 × AC004816.2 expression level) + (1.013 × AC000067.1 expression level) + (0.361 × LINC01235 expression level) + (−2.069 × LINC01010 expression level) + (−0.685 × AL133467.1 expression level) + (0.387 × AC092919.1 expression level) + (−0.676 × AL591468.1 expression level) + (−0.377 × MIR200CHG expression level). The correlation between the 10 prognostic lncRNAs and NETs-related genes was shown in Figure 3F.

Initially, The BC patients in the training and validation cohorts were categorized into high- and low-risk groups according to the training cohort’s medium risk score. The validation cohort was used for internal validation to weight the predictive capability of the model. As represented by the Kaplan–Meier curves, a significantly inferior OS was observed for those in the high-risk group in comparison with those in the low-risk group in both the training (p < 0.001) (Figure 4A) and validation cohorts (p = 0.009) (Figure 4D). The visualized heatmap revealed differential expression of lncRNAs in high- and low-risk groups in training and validation cohorts (Figures 4B, E). Along with the increase in the risk score, the proportion of patients in the high-risk group also increased, and so did the level of mortality (Figures 4C, F).

To further evaluate the possibility that the risk score could serve as an independent BC prognostic signature, the univariate Cox regression analysis was performed by matching the risk score of BC patients in the training cohort with conventional clinicopathological parameters (age, stage, and subtype). In the group of 749 individuals included in this investigation, age, stage, and risk score were associated with the prognosis of the BC patients (Figure 5A). Considering the impact of molecular subtype on BC survival, a multivariable Cox algorithm was conducted to further screen out four prognostic factors of the BC patients (age, stage, subtype, and risk score) (Figure 5B), which was consistent with the result of the univariate Cox regression analysis. The ROC curves and C-index also revealed that the risk score acted as an important role in predicting BC prognosis in training cohort (Figures 5C, D). Based on the four parameters, a nomogram was conducted to predict an individual’s prognosis at 1-, 3-, and 5- years (Figure 5E). The calibration curves showed high consistency between the predicted and the actual 1-, 3- and 5-years survival probabilities (Figure 5F). ROC curves analysis also showed satisfactory AUC values at 1-, 3- and 5- years (0.779, 0.715, and 0.700, respectively) (Figure 5G). In addition, the predictive power and performance of this nomogram was also confirmed in the validation cohort, with better calibration curves and AUC values (Figures 5H, I).

To gain further insight into the specific molecular differences between high- and low-risk groups, differentially expressed genes (as shown in Supplementary Table S4) were identified and functional annotation was conducted by the GSEA. The GSEA with | ES scores | > 0.5, FDR <0.25 and P adjusted value <0.05 for the five most significant pathways showed enrichment of tumor hallmarks in high- and low-risk subgroups. The high-risk group was enriched in HALLMARK_PROTEIN_SECRETION, HALLMARK_GLYCOLYSIS, HALLMARK_G2M_CHECKPOINT, HALLMARK_MTORC1_SIGNALING, HALLMARK_E2F_TARGETS, while the low-risk group was enriched in HALLMARK_ALLOGRAFT_REJECTION, HALLMARK_IL2_STAT5_SIGNALING, HALLMARK_TNFA_SIGNALING_VIA_NFKB, HALLMARK_INTERFERON_GAMMA_RESPONSE, HALLMARK_IL6_JAK_STAT3_SIGNALING (Figure 6A). Then, the potential functions and pathways were identified using GO and KEGG enrichment analyses. GO analysis showed that the differentially expressed genes between the risk score subgroups were enriched in immune response−activating cell surface receptor signaling pathway and immune response−activating signal transduction (Biological Process), immunoglobulin complex and external side of plasma membrane (Cellular Components), antigen binding and immunoglobulin receptor binding (Molecular Function; Figure 6B). KEGG analysis showed that the differentially expressed genes were enriched in cytokine−cytokine receptor interaction, hematopoietic cell lineage, T cell receptor signaling pathway, natural killer cell mediated cytotoxicity, chemokine signaling pathway (Figure 6C). These results indicated that enrichment was mainly focused on immune-related functions and pathways.

The TMB has been implicated as a biomarker for responses to immune checkpoint inhibitor therapy (Chan et al., 2019). The TMB of the high- and low-risk groups of BC patients was shown in Figures 7A, B. The difference in TMB between the two groups was statistically significant (p = 0.022, Figure 7C). Moreover, the patients in the low-TMB group experienced significantly better survival than those in the high-TMB group (p = 0.028, Figure 7D). Moreover, special subgroup analyses stratifying samples were performed according to the combination of TMB status and risk groups. The result showed that some patients with high-risk score in the low- or high-TMB groups had significantly shorter OS than those with low-risk score in the low- or high-TMB groups (p < 0.001, Figure 7E), but no significance could be calculated in the different TMB groups with the same risk scores. Thus, a high level of TMB, commonly considered as a “universal marker,” might fail to accurately predict the reactivity of checkpoint inhibitors across all cancer types. The exploration of the tumor microenvironment (TME) based on the risk predictive model was further conducted, and the proportion of 22 types of tumor-infiltrating immune cells in the high-risk and low-risk groups are shown in Figures 8A, B.

Then, we analyzed the potential correlation between NETs-related lncRNAs and 16 immune cells and the scores of 13 immunological functions using the ssGSEA algorithm. The results revealed that activated dendritic cells (aDCs), B cells, CD8 T cells, dendritic cells (DCs), immature Dendritic Cells (iDCs), mast cells, neutrophils, natural killer cells, plasmacytoid dendritic cells (pDCs), T follicular helper cell, type1 T helper cells, type2 T helper cells, tumor-infiltrating lymphocyte were more predominant in the low-risk group, while regulatory cells (Tregs) were more abundant in the high-risk group (Figure 9A). In addition, immune-related functions were highly enriched in the low-risk group (p < 0.01, Figure 9B). It has also been reported that blocking the immune checkpoint pathway is an extremely promising approach to achieve anti-cancer immunity (Li et al., 2019). Therefore, we made a comparison in the expression discrepancies of checkpoint genes between the high- and low-risk groups. The result was shown in Figure 9C, which indicated that there was significant difference in the expression of all checkpoint genes between the two groups, with all immune checkpoint molecules except CD276 being highly expressed in the high-risk group. These results demonstrate that NETs-related lncRNAs signature can be applied to evaluate the tumor immune microenvironment and the expression of immune checkpoint genes in BC patients.

To assess the prediction performance of NETs-related lncRNAs on drug therapy for BC, we performed an analysis of the relationship between high- and low-risk groups and the efficacy of commonly used therapeutic agents using the “pRRophetic” R package. Our study suggested that the high-risk group showed significant correlation with a higher IC₅₀ for chemotherapeutic and targeted agents such as cisplatin, gemcitabine, paclitaxel, vinorelbine, and gefitinib, which might indicate low-risk group was more suitable for these agents (all p < 0.001, Figures 10A–E). Conversely, lapatinib is more suitable for the high-risk group with a lower IC₅₀ (p < 0.001, Figure 10F). In immunotherapy, TIDE scores were applied to assess the efficacy of immune checkpoint (PD-1 and CTLA-4) inhibitors in the high- and low-risk groups. High-risk patients had significantly lower TIDE scores compared to low-risk patients (p < 0.001, Figure 10G), indicating that patients in the high-risk group might have better response when receiving immune therapy. In the era of immunotherapy, the focus on treatment efficacy should be accompanied by a focus on the management of immune-related adverse effects. The commonly used immunosuppressive agents, such as methotrexate and rapamycin, differed significantly between the two groups, with the high-risk group having a higher IC₅₀ (p < 0.001, Figures 10H, I).

First, the paired lncRNAs expression profiles of breast samples were analyzed by TCGA databases. Compared with adjacent normal breast tissues, the lncRNAs expression of SFTA1P, ACTA2-AS1, AC004816.2, LINC01235, and AL133467.1 were lower, while the expression of AC000067.1, AL591468.1, and MIR200CHG was higher in BC tissues (Figure 11A). Only LINC01010 and AC092919.1 were not statistically significant in paired tissues. To further confirm the validity of this model, the expression levels of all lncRNAs in 10 BC tumor tissues and matched normal tissues from our hospital were detected using qRT-PCR. Except for LINC01010, their expression trend was consistent with the bioinformatic analysis results (Figure 11B). The qRT-PCR results indicated that our bioinformatics analysis was accurate and reliable, reinforcing our conclusions from these data.