Breast cancer gene expression data, including coding and non-coding RNA sequence data, were acquired from TCGA (https://cancergenome.nih.gov/) together with the corresponding clinical information. Until 2017, 1,098 breast cancer samples were available from TCGA, though in the present study only those including patient survival status were selected (n=768); this enabled the determination of any association between the expression of lncRNAs of the lncRNA-expression signature and the corresponding OS time for breast cancer. These 768 breast cancer samples were divided equally into a training set (to identify the gene expression signature) and a validation set (to validate the gene expression signature). To confirm the expression levels of the differentially expressed genes, the gene expression dataset GSE5764 (21), containing 10 breast cancer tissue samples and 20 non-cancerous samples, was downloaded from GEO (https://www.ncbi.nlm.nih.gov/geo/; Affymetrix GPL570 platform, Affymetrix Human Genome U133 Plus 2.0 Array; Affymetrix; Thermo Fisher Scientific, Inc.).

RNA genes downloaded from TCGA were compared with published lncRNAs from the MiTranscriptome database (http://mitranscriptome.org/). Potential lncRNAs were identified as transcriptome sequences that were mapped to corresponding lncRNAs, rather than any protein-coding region, and were not identified as protein-coding genes in the National Center for Biotechnology Information database (https://www.ncbi.nlm.nih.gov/).

Raw read counts of the transcriptomic data from TCGA were normalized using the quartile normalization method and logarithmically transformed to a normal distribution. The Bioconductor package DESeq2 (https://www.bioconductor.org/packages/release/bioc/html/DESeq2.html, version 1.24.0) was used to perform the normalization and identify differentially expressed genes in breast cancer samples compared with adjacent normal tissues, with an adjusted P<0.05 and an absolute log2-based fold-change >0.5. For gene expression data from GEO, the R package limma (https://www.bioconductor.org/packages/release/bioc/html/limma.html, version 3.40.6) was used to conduct differential expression analysis.

To establish a prognostic signature for breast cancer, a two-step method was employed using the R package SIS (https://CRAN.R-project.org/package=SIS, version 0.8–6) for sure independence screening. Firstly, univariate Cox regression analysis was performed to identify survival-associated genes. Secondly, SIS (based on the least absolute shrinkage and selection operator, Cox-penalized regression model) was used to identify important variables and construct multi-gene-based prognostic signatures for OS rate prediction.

To identify genes that correlated with the four lncRNAs of the prognostic signature, data from TCGA were used to evaluate the pairwise Pearson's correlation between the expression levels of the target lncRNAs. Only associated genes with an absolute r≥0.3 and a significant correlation (P<0.05) were retained. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG; http://www.genome.jp/kegg/) pathway analyses were performed using the Database for Annotation, Visualization and Integrated Discovery (DAVID; http://david.ncifcrf.gov/).

Kaplan-Meier analysis was used to estimate the performance of the prognostic signatures, and log-rank test was performed to evaluate statistical significance. A risk score was calculated for each patient according to the formula of the four-lncRNA signature, and patients were divided into high- and low-risk groups using the median score as a cut-off. Receiver operating characteristic (ROC) analysis was used to evaluate the sensitivity and specificity of the four-lncRNA signature and other biomarkers, including TP53, MKI67, ESR1, PGR, ERBB2 and HOTAIR. All statistical analyses were conducted using R 3.5.2 (https://www.r-project.org/), and P<0.05 was considered to indicate a statistically significant difference.

All 768 patients were diagnosed with breast cancer based on clinicopathological evaluation. The clinical stage and histological subtype were determined using the Tumor-Node-Metastasis staging (22) and immunohistochemical molecular typing methods, respectively. According to the data, the estrogen receptor (ER), progesterone receptor (PR) and human epidermal growth factor receptor 2 (HER2) status of each patient was indicated as positive, negative or indeterminate. In addition, the ranges of the OS and relapse-free survival times were 1–8,605 and 1–8,556 days, respectively. Patient characteristics are displayed in Table I.

Differential expression analysis was employed to select differentially expressed RNAs in normal and cancerous tissues; a total of 8,854 upregulated and 5,939 downregulated RNAs were identified. Subsequently, all possible combinations of four lncRNAs were analyzed and compared using the two-step Cox regression method. A total of 7 models consisting of four lncRNAs were identified (Table SI); among these candidates, one model was identified as the most suitable for predicting the OS of patients with breast cancer. Patients were divided into high- and low-risk groups using the median risk score as a cut-off, with the risk score calculated as follows: Risk score=(−0.015× expression value of PVT1) + (−0.193× expression value of MAPT-AS1) + (−0.116× expression value of LINC00667) + (0.098× expression value of LINC00938). The coefficients of this formula were derived from multivariate Cox regression analysis (Table SI). Kaplan-Meier analysis was also performed to determine the association between the expression levels of the four-lncRNA signature and patient OS. Compared with those of the low-risk group, high-risk patients exhibited significantly poorer OS rate (log-rank P=0.0151; Fig. 1A).

To confirm the predictive capacity of the four-lncRNA signature identified in the training set, the equivalent analyses were also performed in the validation set. Patients were divided into low- and high-risk groups, and the differences between patient OS rates were compared using Kaplan-Meier analysis. Patients in the high-risk group possessed significantly lower OS rate than those of patients in the low-risk group (log-rank P=0.0016; Fig. 1B), which was consistent with the findings from the training set. Furthermore, ROC analysis was performed to evaluate the sensitivity and specificity of survival prediction; the area under the curve (AUC) was 0.641 (Fig. 2A), indicating that the four-lncRNA signature was able to accurately predict the survival of patients with breast cancer.

ROC analyses were performed to investigate whether the four-lncRNA signature was applicable to different breast cancer stages and molecular subtypes. In stages I–IV, the AUC values were 0.595, 0.687, 0.634 and 0.645, respectively (Fig. S1), indicating that the four-lncRNA signature was able to predict the survival of patients at different clinical stages of breast cancer. Regarding subtype, the AUC values in the four molecular subgroups were 0.637, 0.654, 0.688 and 0.613, respectively (Fig. S2), suggesting that the four-lncRNA signature served as a prognostic indicator for patients with different breast cancer subtypes.

For further clarification, the performance of the four-lncRNA signature was compared with that of several known breast cancer biomarkers, including TP53, MKI67, ESR1, PGR, ERBB2 and HOTAIR, using ROC and Kaplan-Meier analyses. The sensitivity and specificity of these six known biomarkers are displayed in Fig. 2B. The AUC values of TP53, MKI67, ESR1, PGR, ERBB2 and HOTAIR were 0.574, 0.483, 0.510, 0.501, 0.501 and 0.473 (data not shown), respectively, while the AUC value (0.641) of the four-lncRNA signature was greater. Kaplan-Meier analysis revealed that only TP53 was significantly associated with patient OS (Fig. 3A; log-rank P=0.0396), while the other selected biomarkers were not (Fig. 3B-F). Moreover, the four lncRNAs from the identified model were also evaluated. ROC curve analysis generated AUC values for PVT1, MAPT-AS1, LINC00667 and LINC00938 as 0.532, 0.553, 0.550 and 0.480, respectively (Fig. S3), which were smaller than that of the four-lncRNA signature. In addition, the results of Kaplan-Meier analysis indicated that the differential expression of these lncRNAs was not significantly associated with the OS rate of patients with breast cancer (log-rank P>0.05; Fig. S4).

To further investigate the potential functions of the four lncRNAs of the signature, gene expression data from the GSE5764 dataset were downloaded from the GEO database, and differential expression analysis was performed. The fold-change values of PVT1, MAPT-AS1, LINC00667 and LINC00938 were 2.031, 3.057, 1.579, 0.455, respectively. This result indicated that PVT1, MAPT-AS1 and LINC00667 were upregulated in breast cancer tissues, and that LINC00938 was downregulated. GO enrichment and KEGG pathway analyses were conducted. According to the results of GO analysis, the primary PVT1-associated functions were ‘transcription elongation’, ‘mitochondrial electron transport’ and ‘endoplasmic reticulum-associated degradation’ (Fig. 4A). MAPT-AS1 was associated with ‘cilium morphogenesis’ and ‘cilium assembly’, ‘intraciliary retrograde transport’ and ‘neurological system process’ (Fig. 4B). LINC00667 was associated with ‘regulation of transcription, DNA-templated’, ‘centrosome organization’ and ‘regulation of cell morphogenesis’ (Fig. 4C), and LINC00983 was associated with biological processes involved in ‘cilium assembly’ and ‘S-adenosylmethionine metabolic process’ (Fig. 4D). Moreover, KEGG analysis also indicated several biological processes and pathways that potentially associated with the signature lncRNAs, including ‘Notch signaling pathway’, ‘oxidative phosphorylation’ and ‘Huntington's disease’ (Fig. 4E).