Novel lncRNAs with diagnostic or prognostic value screened out from breast cancer via bioinformatics analyses

Transcriptome sequencing data with clinical information collected from databases of breast cancer tissues

Using the TCGAbiolinks (Colaprico et al., 2016) R package, the breast cancer subdata of the transcriptome sequence from TCGA (the HTSeq-Counts type) were downloaded, resulting in a total of 1,222 breast tissue samples. The data from 794 cases of diagnosed invasive ductal carcinoma in breast tissues and 92 cases with normal tissue were also collected and downloaded with their clinical information.

Screening of differentially expressed lncRNAs and mRNAs in breast cancer

The Deseq2 of R package was used to screen out the differentially expressed genes (DEGs) between the breast cancer and paracancerous tissue groups. |Log2 fold change| ≥ 1 and adjust p value <0.05 were applied as the selected threshold, and then the annotations were performed based on those candidates. The annotation films of Gencodev27.gtf in the GENCODE database (https://www.gencodegenes.org) served to distinguish the lncRNAs from RNA molecules.

lncRNA sequencing data with correlation analyses between extracellular vesicles in plasma and tissues of breast cancer patients

The plasma samples of breast cancer and benign breast tumor groups (20 cases/group) for EV extracts were collected and treated from the 6th Affiliated Hospital of Shenzhen University, according to documents approved by the Medical Ethics Committee of the Union Shenzhen Hospital, Huazhong University of Science and Technology (0720001, ky-2020-039-02), and clinical trial registration (ChiCTR19002505). A written informed consent was obtained from every participant for this study.

Before RNA-seq, 10 plasma samples were pooled into one with an equal volume from same group, and the total RNAs of plasma were extracted using the exoRNeasy Serum/Plasma Maxi Kit (QIAGEN, USA). The large RNA molecules were sequenced to detect and analyze the expression levels of the lncRNAs.

The data were also compared with the expression of plasma EVs in other cancers from the exoRbase in order to find the disease-related lncRNAs in breast cancer. The lncRNA candidates of plasma EV were compared with lncRNAs of tissues, in order to screen out cancer-related lncRNAs, which were both differentially expressed in plasma EVs and breast tissues.

Construction of weighted gene co-expression network

A weighted gene co-expression network analysis (WGCNA) was conducted (Langfelder & Horvath, 2008) to investigate the function of the lncRNAs by extracting their expression data from the TCGA and obtaining the clinical information from the samples. The data were processed and the function “goodSamplesGenes” was used to test the gene expression matrix to confirm whether there was a missing value. The “hclust” function was used to cluster the samples to remove the outlier samples. Using the Pearson correlation coefficient, the linear correlation degree between the two genes was analysed and the gene expression matrix was transformed into the gene relationship matrix. Then the inter-gene relationship matrix was transformed into a weighted scale-free network; the index β (soft threshold) was taken for the previously obtained Pearson correlation coefficient between genes, and the β value was adjusted to make the network meet the scale-free characteristics, resulting in the adjacency matrix of genes. The adjacency matrix was upgraded to a more rigorous and reliable topological overlap matrix (TOM). Using the difference degree of TOM as the clustering distance for cluster analysis, the last branch of the cluster tree represented genes, and the genes with a high overlap were clustered together to form large branches. Using the dynamic tree pruning (DynamicTreeCut), all genes were divided into different modules, and the number of genes contained in the module were not less than 30, and the modules were distinguished by colour (Langfelder, Zhang & Horvath, 2008). The genes in the module were highly related, and there were a small number of genes, which reflected the expression characteristics of the whole module. This was called the feature vector gene (Module eigengene, ME) of the module.

Identification of the key lncRNA with clinical information

After the module identification was completed, it was combined with clinical information for analysis. The operation of WGCNA was based on continuous variables and classification variables, such as sample type (tumour or normal). The TNM was set to continuous variables. The sample type was converted to 0 or 1 min and the TNM staging was converted to the corresponding number (0, 1, 2, 3, 4, 5). If there was no clear clinical information in the sample, the samples were treated as a null value. Other clinical data such as age and survival time were not modified. The correlation between the modules and clinical features was achieved by calculating the Pearson correlation coefficient between ME and clinical features. Additional studies were performed to analyse the gene composition in the selected modules, identify the lncRNAs and the top five mRNAs related to lncRNAs, use Cytoscape to visualize the results, and apply g:profiler (Raudvere et al., 2019) for enrichment analysis to predict the potential biological functions of lncRNA in breast cancer.

Assessment of biomarker function of core lncRNAs

The ROC curve was used to evaluate the diagnostic ability of the key lncRNAs to distinguish between breast cancer samples (n = 794) and normal samples (n = 92); the standard of screening was AUC area >70%. Univariate cox regression analysis was used in conjunction with the clinical information to study the relationship between the initial age at diagnosis, tumour stage, expression of key lncRNA, and total survival time. The initial age at diagnosis was categorized into a young and old group according to the median age. According to the tumour stage, stage I and II were classified as early stage, and stage III and greater were classified as late stage. The selected factors were then analysed by multivariate cox regression analysis to screen the factors with independent prognostic function.

Partial verification of key lncRNA

The expression data of EV RNA (118 cases of normal and 140 cases of tumour) and tissue RNA (12 cases of normal and 74 cases of tumour) of breast cancer were downloaded from exoRBase and Gene Expression Omnibus database (GEO). Log ₂ (TPM+1) was used as the expression of key lncRNA, and the expression of key lncRNA was analysed using the Student’s t-test. ROC curves were then used to analyse the diagnostic roles of key lncRNAs in EVs and tissues, respectively.

Statistical analysis

All data statistics were carried out by R software (version 3.5). The screening of differentially expressed genes (DEGs) was completed by R packet Deseq2. Genes with an |Log2 fold change | ≥1 and adjusted p-value <0.05 were considered to be differentially expressed. The resulting adjusted p-value was obtained using the Benjamini and Hochberg’s approach. The diagnostic ability and prognostic ability of lncRNA were performed by ROC curve (R packet pROC) and survival analysis (R packet survival and survminer), respectively. Univariate cox regression analysis and multivariate cox regression analysis were used to find lncRNA with independent prognosis; HR ≠ 1 and p < 0.05 were used as screening factors. In the verification analysis, the expression of lncRNA was treated with log2 (TPM+1), and the comparison of expression levels between normal tissues and tumors was analyzed using the Student’s t-test.

Differential expression and gene annotation of lncRNAs in breast cancer tissues

A total of 11,321 differentially expressed genes (|log2Fold Change| ≥ 1, adjusted p-value <0.05) were obtained from 56,512 genes of 886 tissue samples using DESeq2. These were combined with the annotation information in GENCODEv27.gtf and 4,879 differentially expressed lncRNAs and 5,304 differentially expressed mRNAs were obtained. Further study showed that 2,883 of these differentially expressed lncRNAs were up-regulated and 1,996 were down regulated. Of the differentially expressed mRNAs, 3,226 were up-regulated and 2,078 were down-regulated. The different expression patterns of lncRNAs and mRNAs are shown in Fig. 1.

Comparative analysis of plasma lncRNAs in plasma EV of breast cancer patients

According to the threshold value of | log₂Fold Change | ≥1 and adjust p value <0.05, a total of 2,714 DEGs were screened using DESeq2 software. Then 155 lncRNAs and 2,450 mRNAs were selected by annotating the differentially expressed genes. Among the differentially expressed lncRNAs, 72 were up-regulated and 83 were down-regulated. In addition, among the differentially expressed mRNAs, 1,619 were up-regulated and 831 were down regulated (Fig. 2). The top 20 differentially expressed lncRNAs are shown in Table 1. A total of 12 samples of blood extracellular vesicles (GSE100063) from patients with colon cancer (CRC), 21 samples of blood extracellular vesicles (GSE100207) from patients with hepatocellular carcinoma (HCC), and 12 samples of blood extracellular vesicles (GSE100207) from patients with pancreatic cancer (PAAD) were downloaded from the exoRBase database. Furthermore, the expression profiles of these samples were analysed using Limma package software, and lncRNAs were re-annotated by GENCODEv27.gtf. Thirty-seven differentially expressed lncRNAs were found in CRC, of which 33 were up-regulated and four were down-regulated; sixty-nine were found in HCC, of which 21 were up-regulated and 48 were down-regulated; a total of 28 were found in PAAD, of which 25 were up-regulated and three were down-regulated. A comparison of lncRNAs related to these diseases with those of the breast cancer lncRNAs showed that 149 lncRNAs were highly correlated with breast cancer. The data were further compared with lncRNAs in tissue expression profiles and 41 lncRNAs were found to be differentially expressed in both breast cancer tissues and plasma EVs.

Construction of weighted correlation network

We extracted information on the expression of 41 lncRNAs and clinical information from the TCGA database for 866 tissue samples. This was integrated with the data previously obtained for 5,304 mRNA from the TCGA tissue samples to construct a weighted correlation network. The screening range threshold was defined as 0 to 30 and the correlation degree of log (k) and log (p (k)) under each threshold, as well as the corresponding average connectivity and average correlation degree were calculated. The soft threshold β = 3 met the requirements of scale-free network (Fig. S1). Therefore, the soft threshold = 3 was selected to calculate the adjacency matrix. The adjacency matrix was promoted to the TOM matrix, and the degree of TOM dissimilarity was used as the clustering distance for clustering analysis. A total of 5,345 genes were divided into different modules. As the number of genes contained in the module were not less than 30, and these modules were distinguished by the colour, a total of 19 gene modules were obtained. The Pearson correlation coefficient analysis module was used to analyse the correlation degree of each clinical feature and disease. Five modules (blue–green module, yellow module, purple module, blue module, and brown module) were found to have a relatively high correlation with the T stage and the TNM comprehensive stage of breast cancer (Fig. 3).

Composition and function analyses of important modules

The blue–green, yellow, purple, blue and brown modules related to TNM staging closely were analysed, and the lncRNAs and mRNAs with similar expression in the same module were identified, with a total of 28 of lncRNAs and 3,901 of mRNAs (Table 2). g: profiler was used to analyse the functional enrichment of important modules to explore their important role in breast cancer (Fig. 4). Among them, the genes of blue–green module were involved in the interaction of neuroactive ligand receptors and played roles in the cAMP signalling pathway, transcriptional imbalance in cancer, peroxisome proliferator activated receptors (PPARs) pathway, and PI3K Akt signalling pathway. The genes of yellow module might be involved in the formation of adhesive plaque and the interaction of extracellular matrix (ECM) receptors; the genes of purple module gene might be related DNA replication, cell cycle regulation, homologous recombination and the p53 signalling pathway related to the occurrence and development of breast cancer; while the genes of brown module might be involved in pyrimidine nucleotide metabolism and other processes. In the blue–green module, 19 lncRNAs and 2,273 mRNAs were obtained; in the yellow module, two lncRNAs and 289 mRNAs were obtained; in the purple module, 108 mRNAs were obtained; in the blue module, four lncRNAs and 658 mRNAs were obtained; in the brown module, three lncRNAs and 573 mRNAs were obtained.

Furthermore, the candidate lncRNAs in the modules were further analysed and sorted by TOM value to find the mRNAs with a relatively high TOM value, which were the potential targeted regulatory genes of the lncRNAs. The top five were selected (Table 3), and the co-expression network was constructed according to the regulatory relationship among of them (Fig. 5). In the co-expression network, TGFBR2, CAV1, PDE2A, LDB2, EBF1, and other key genes were regulated by multiple lncRNAs.

Capability analysis of key lncRNAs as biomarkers

In order to explore the diagnostic ability of the 28 candidate lncRNAs in breast cancer, the ROC curve was analysed. An AUC value >70% was set as the screening criteria; a total of eight lncRNAs were found to be able to be used for diagnosing breast cancer. Among them, the AUC values of LINC00514 (AUC = 86.7928%), C15orf54 (AUC = 86.8365%), WFDC21P (AUC = 77.7562%), AL157935.1 (AUC = 76.484%), AC124798.1 (AUC = 75.0541%), AL136162.1 (AUC = 70.6103%), LINC01117 (AUC = 70.9059%), SNHG3 (AUC = 70.6789%) were >70% (Fig. 6). In combination with the clinical information, the age of initial diagnosis, tumour stage, and the role of AL355974.2 (HR = 0.79, p = 0.0077) in prognosis were determined using univariate cox regression analysis (Table 4). These factors were analysed by multivariate cox regression analysis, and it was found that the early diagnosis of the tumour (young diagnosis age, early tumour stage) was of great significance for improved survival. AL355974.2 (HR = 0.8103, p = 0.02298) could also be used as an independent prognostic factor, and as a protective factor as its high expression helps to maintain the survival rate of patients. Our results were statistically significant (HR ≠ 1 and p <0.05), as shown in Fig. 7.

Extending verification of key lncRNAs

In order to further confirm the function and role of these lncRNA, we downloaded the RNA expression data in EVs and tissues of breast cancer from exoRBase (GSE93078) and GEO (GSE134359). In the EV dataset, seven lncRNAs expression data, including LINC00514, C15orf54, AL157935.1, AC124798.1, AL136162.1, LINC01117 and SNHG3 were extracted (Fig. 8). The expression levels of C15orf54, LINC01117 and SNHG3 varied between normal and tumour groups, and LINC01117 showed an up-regulated trend in tumours, while C15orf54 and SNHG3 showed a down-regulated trend in tumours. The AUC value of the seven lncRNAs were all less than 70% (Table 5), indicating that these lncRNAs do not have good diagnostic ability in the EV for breast cancer.

In the tissue data set, the expression data of five lncRNAs, such as LINC00514, C15orf54, AL157935.1, LINC01117 and SNHG3 were extracted, and the expression levels are shown in Fig. 9. Among them, the expression levels of C15orf54, AL157935.1, LINC01117 and SNHG3 were different between normal and tumour groups, and all of them were up-regulated in tumours. The ROC curves of five lncRNAs were drawn, as shown in Fig. 10. LncRNAs such as C15orf54 (AUC = 100%), AL157935.1 (AUC = 99.4369%), LINC01117 (AUC = 76.3514%) and SNHG3 (AUC = 88.4009%) had good diagnostic ability in tissue samples.