In order to investigate the role of aberrant expression of lncRNAs in breast cancer, a bioinformatics approach was used to predict the differentially-expressed genes. The symbols and names of the lncRNA genes were obtained from the database of the HUGO Gene Nomenclature Committee (genenames.org). The breast invasive carcinoma dataset was obtained from The Cancer Genome Atlas (TCGA) at the cBioPortal for Cancer Genomics, consisting of 1,105 samples (cbioportal.org). By selecting the cancer study and genomic profiles, the genes of lncRNAs that have been downloaded may be entered and the information submitted. miRNAs may be selected in the same way.

In order to identify the miRNAs which are able to target lncRNAs, the binding of lncRNAs to miRNAs (including the folded RNA predicted structure of the lncRNAs and miRNAs) was predicted using the bioinformatics tool RegRNA 2.0 (regrna2.mbc.nctu.edu.tw/detection.html). The protein sequence was obtained from the GenBank database of the National Center for Biotechnology Information (NCBI; ncbi.nlm.nih.gov/genbank). In addition, when predicting the miRNA target sites, the minimum folding free energy was set under <-25 and the system score was set to >160. An increased score indicates a stronger binding ability. Due to the large number of lncRNAs which are differentially expressed in breast cancer, the lncRNAs selected fell above 15% alteration frequency. Additionally, the lncRNA sequences associated with Homo sapiens were searched in advance using the NCBI database.

In order to predict which genes to target with screening miRNAs, three databases were used: TargetScan (targetscan.org); Microcosm Targets (ebi.ac.uk/Enright-srv/microcosm/htdocs/targets/v5); and PicTar (pictar.mdc-berlin.de/cgi-bin/PicTar_vertebrate.cgi). Genes which were identified using all of the three databases were used, in order to limit the number of false positive results. In addition, it was confirmed that the lncRNAs exhibited alterations >15%, in order to maximize the clarity of the network diagram. The lncRNA-miRNA-mRNA interaction network was constructed using Cytoscape software (version 3.0; http://www.cytoscape.org/download.php).

In order to investigate the association between clinical pathological features and lncRNAs, computer analysis based on the RPISeq RNA-protein interaction prediction program (pridb.gdcb.iastate.edu/RPISeq), which uses random forest (RF) and support vector machine (SVM) classifiers to predict the association between lncRNAs and breast cancer-associated proteins. The protein sequences were obtained from the NCBI GenBank database.

In order to further investigate the biological effects of aberrantly-expressed lncRNAs and miRNAs in breast cancer, GO enrichment of the target gene were carried out using the GOrilla tool (cbl-gorilla.cs.technion.ac.il). For each GO term, a list of associated genes is returned with the most optimal at the top of the list. Each gene name is specified by the gene symbol and followed by a short description of the gene.

In the present study, among 2,772 lncRNAs analyzed, 626 lncRNAs exhibited alterations in 1–22% of breast cancer cases, with the exception of 22 invalid gene symbols. A total of 41 lncRNAs exhibited alterations in ≥10% cases, and 15 of them (CASC8, PVT1, CCAT1, CCDC26, LINC00536, LINC00861, LINC00964, PCAT1, PCAT2, LINC0097, HPYR1, BAALC-AS2, LINC00051, LINC00535) accounted for a total of 30%. The 15 identified lncRNAs with alterations in breast cancer have been well-studied previously. Zhang et al (12) evaluated the expression of CCAT1 by reverse-transcription quantitative polymerase chain reaction (RT-qPCR) analysis, demonstrating that the expression of the lncRNA CCAT1 is increased in breast cancer tissue compared with healthy adjacent tissues, and that it may be a novel biomarker of poor prognosis in patients with breast cancer. For the miRNAs, it was identified in the present study that the majority of them exhibited alterations in the 971 records of breast cancer data; specifically, 1,172/1,189 miRNAs exhibited alterations and the frequency ranged between 1 and 21%. miR-510 overexpression in breast cancer has been demonstrated to increase cell growth, migration and invasion (13). The results of the present study are presented in Fig. 1.

Although it is known that miRNAs are able to target a number of protein-coding genes, little is known about whether miRNAs are able to target lncRNAs. In order to elucidate the precise mechanism underlying the role of lncRNAs and miRNAs in breast carcinoma, RegRNA 2.0 was used to analyze the functional association between lncRNAs and miRNAs. RegRNA 2.0 is an integrated web server used to compare the homologs of regulatory RNA motifs and elements against an input mRNA sequence. The predicted miRNAs from this software ought to intersect with the breast invasive carcinoma dataset from TCGA. Among the significantly differentially-expressed lncRNAs, using a threshold alteration frequency >10%, 393 miRNAs are able to exert regulatory functions on 41 lncRNAs. The lncRNAs with the highest alteration frequencies were selected, demonstrating that 13 lncRNAs, targeted by 158 miRNAs, are associated with breast cancer. It was predicted, for example, that lncRNA CCAT1 may be regulated by hsa-miR-595 and hsa-miR-345-5p, and linc00536 may be regulated by hsa-miR-93-5p and hsa-miR-214-5p. By contrast, hsa-miR-93-5p may be targeted by linc00535 and linc00964. A previous study demonstrated that miR-21 is able to regulate lncRNA expression. RT-qPCR analysis of lncRNA demonstrated that miR-21 is capable of downregulating the lncRNA GAS5 (14).

The folded RNA predicted structure of the lncRNAs and miRNAs was analyzed using RegRNA2.0 software, and the partial reliability information of pair probabilities are depicted in Fig. 2. Hairpin loops are the most frequently-predicted secondary structure, which may facilitate binding of lncRNAs to miRNAs. lncRNA CCAT1 is able to bind to hsa-miR-1178 and hsa-miR-1304-3p, and linc00861 is able to target hsa-miR-510. Additionally, previous studies demonstrated that lncRNAs harbored putative miRNA regulatory elements and served important roles in the miRNA regulatory network. Pilyugin and Irminger-Finger (15) observed that the newly-discovered lncRNA BARD1 9′L shared miR-203 and miR-101 miRNA response elements (MREs) with BARD1 mRNA in their homologous 3′ untranslated regions. The network constructed in the present study demonstrated the unknown complexity of non-coding RNA regulatory interactions and how lncRNAs may serve as important factors in the miRNA regulatory network, with the aim of elucidating the role of these interactions in disease processes.

A previous study demonstrated that miRNAs served a role in regulating gene expression by causing the degradation of target mRNAs. Lin et al (16) demonstrated that miR-33b inhibits breast cancer metastasis by targeting high mobility group AT-hook 2, spalt like transcription factor 4 and twist family bHLH transcription factor 1. In the present study, 35 miRNAs were identified to be differentially-expressed in breast cancer. In order to eliminate false positive rates of the target prediction, only the miRNA-mRNA pairs simultaneously predicted by ≥2 applications were taken forward. A total of 549 genes were predicted to be targets of aberrantly expressed miRNAs. The majority of the targets were associated with cancer, including cell division cycle 34 (CDC34), mitogen-activated protein kinase kinase kinase 11 (MAP3K11), AT-rich interaction domain 4B and caspase recruitment domain family member 11. The target genes exhibited functions in cell proliferation, apoptosis, invasion and metastasis. The genes were submitted to Cytoscape for visualization (Fig. 3).

A miRNA-mRNA regulatory network was constructed in the present study. As presented in the network diagram (Fig. 3), a gene may be targeted by multiple miRNAs. The number of target genes of hsa-let-7e-5p is 10 (solute carrier family 35 member D2, CDC34, basic leucine zipper and W2 domains 1, adrenoceptor β-2, heart and neural crest derivatives expressed 1, glutaminyl-tRNA synthetase, adaptor related protein complex 1 sigma 1 subunit, G protein subunit γ-5, leucine rich repeats and immunoglobulin like domains 3 and coilin). Apolipoprotein C2 mRNA is able to be co-regulated by hsa-miR-637, hsa-miR-422a and hsa-miR-765. A previous study demonstrated the interaction between miR155 and tumor protein p53 inducible nuclear protein 1 (TP53INP1). TP53INP1 is a pro-apoptotic stress-induced p53 target gene. Gironella et al (17) demonstrated that TP53INP1 was a target of miR-155 using bioinformatics and cell experiments. In the present study, the regulation of miR155-TP53INP1 was demonstrated in breast cancer. Zhang et al (18) demonstrated that overexpression of miR-155 led to the downregulation of TP53INP1 and reversed the effect of TP53INP1, by promoting breast cell proliferation and suppressing cellular apoptosis. It has additionally been demonstrated that neurofibromatosis type 1 (NF1β) is targeted by miR-106a and miR-21, and that the NF1β mRNA flanking region around the core sequence for miR-21 recognition is conserved in vertebrates up to 500 bp (19). NF1β knockdown led to a phenocopy of overexpressing miR-21 in the study by Dellago et al (20). In the present study, it was hypothesized that certain common target mRNAs may be detected between different dysregulated miRNAs. From the miRNA-mRNA regulatory network constructed in the present study, the regulation of hsa-miR-199a-5p-MAP3K11 was observed. Byrnes et al (21) demonstrated that miR-199a-5p acts as a tumor suppressor in carcinogenesis, and that its downregulation contributes to enhanced cellular proliferation by targeting MAP3K11. Further details of the results of the present study are presented in Fig. 3.

The regulatory network constructed in the present study exhibited 601 nodes and 706 edges, demonstrating the complexity regulatory interactions between lncRNAs, miRNAs and target genes. However, the functional roles of a number of lncRNAs and miRNAs in the regulatory network have been rarely reported in breast cancer. Cesana et al (22) demonstrated that a muscle-specific lncRNA, linc-MD1, ‘sponges’ miR-133 and miR-135 to regulate the expression of mastermind like transcriptional coactivator 1 and myocyte enhancer factor 2C, transcription factors that activate muscle-specific gene expression. The results of this previous study indicated that there may be direct competition for miRNA binding between lncRNAs and mRNAs, and that this serves a role at the post-transcriptional level.

In order to further investigate the regulatory information network associated with the clinical pathological features of breast cancer, RPISeq software was used to predict the potential target proteins of lncRNAs. Interaction probabilities generated by RPISeq range between 0 and 1. In performance evaluation experiments, predictions with probabilities >0.5 were considered ‘positive,’ indicating that the corresponding RNA and protein are likely to interact. The score of the predicted interaction possibility between lncRNAs and clinical pathological features are presented in Fig. 4. The results of the present study demonstrated that a number of lncRNAs, including PVT1, CCAT1 and linc00861, exhibited a strong possibility of interaction with receptor tyrosine protein kinase erbB-2 (HER2), estrogen receptor (ER), progesterone receptor (PR), breast cancer type 1 susceptibility protein, breast cancer type 2 susceptibility protein, cellular tumor antigen p53, phosphatase and tensin homolog, and tumor necrosis factor, with scores of RF>0.6 and SVM>0.9. It was hypothesized that lncRNAs PVT1, CCAT1 and linc00861 may govern fundamental biological processes and may be used as clinical biomarkers for diagnosis. Merry et al (23) reported that HER2 amplification affects the expression of lincRNAs, using RNA sequencing in human breast cancer. The results of the present study may provide a foundation for future experiments, leading to increased understanding of the associations between RNA-based biomarkers and clinical pathological features.

In order to investigate the role of lncRNAs and miRNAs in breast cancer, a systems biology approach was used to examine the functions of predicted targets. Genes (548) were recognized by gene symbol out of 549 gene terms entered. A total of 442 of these genes were associated with a GO term. The results of the enriched biological processes and molecular functions analysis, which were generated with the GOrilla tool, are presented in Fig. 5. GOrilla is a widely-used application for identifying and visualizing enriched GO terms in ranked lists of genes (24). Each differential target gene may return a P-value; if this value is <5%, it denotes an enrichment gene. Enrichment genes may contribute to multiple biological processes, including the type I interferon (IFN) signaling pathway, natural killer cell activation in the immune response, the response to double stranded RNA and B cell proliferation. The most notable biological function was the regulation of the type I IFN signaling pathway. Browne et al (25) demonstrated that type I IFN, including IFNα, inhibit the replication of human immunodeficiency virus I by upregulating the expression of important genes. An additional enrichment ranking of GO and signaling pathways indicated that cytokines may be the most enriched genes. A study of eight cytokine genes conducted by Kim et al (26) demonstrated that 74 (24%) patients were classified using eight cytokine genes as exhibiting prevalent depression, and 19 (8%) and 25 (10%) patients were classified with persistent and incident depression, respectively. These previous results supported the role of cytokines in the etiology of depression associated with breast cancer.

The functional analysis conducted in the present study was supported by the above previous studies, and demonstrated using a variety of bioinformatic approaches, that lncRNAs and miRNAs exerted regulatory effects on breast cancer activation, by affecting molecular functions and signaling pathways.