RNA-sequencing data, updated clinical data, and sample information of breast invasive cancer (BRCA) cohort were downloaded from the TCGA data portal (https://tcga-data.nci.nih.gov/tcga/) and 107 patients in GSE58812 for external validation was obtained from the Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/) database. Patients who fulfilled the following selection criteria: (1) diagnosed with histologically confirmed BRCA; (2) available RNA-seq data; (3) available follow up were eligible for study enrollment. Finally, 1054 patients from TCGA-BRCA, 107 patients from GSE58812 were included for further analysis. For TCGA-BRCA data analysis, the established RNA matrix file was merged with Ref-Seq and used for all subsequent RNA-seq expression analyses. For the data from GEO datasets, the probe ID for each gene was converted to a gene symbols. All raw data were analysed using the ‘limma’ package in R (version 4.2.1, http://www.bioconductor.org/).

In the current study, we established 1542 encoding RNA-binding proteins (RBPs) [9]. Correlations between lncRNAs and RBPs were calculated using Pearson’s correlation coefficient using the limma R package [18]. A lncRNA with a correlation coefficient $|$R${}^2$$|$ $>$0.3 and p 0.001 was considered to be an RBP-related lncRNA. An RBPLSig was established based on a linear combination of the regression coefficient derived from the multivariate Cox regression model coefficients ($\beta$) and expression levels of the RBP-related lncRNAs. To calculated the risk score, the following formula was used: ${\sum }_{i=1}^n\beta i*\left(\text{expression of lncRNA}\right).$Using the median score of the RBPLSig, patients were divided into high- or low- groups. Survival curves analysis was performed between high- or low- groups by using the “survival” and “survminer” packages. Besides, the receiver operating characteristic (ROC) analysis was performed by using the “timeROC” and “survivalROC” R packages to assess the prognostic accuracy of RBPLSig value. The distribution of patients with different risk scores was investigated by principal components analysis (PCA) using the R package “princomp”, and the PCA 3D graph were visualized using the R packages “catterplot3D”. The co-expression network and Sankey diagram were used to established and visualized the correlation between the poor prognostic RBP-related-lncRNAs in the prognostic model and their co-expressed RBP mRNAs.

To better research and further explore the molecular mechanisms underlying the prognostic difference between high- and low-risk groups the biological function of GSEA was performed to by using the GSEA tools (http://www.broadinstitute.org/gsea) [19]. Significant predefined biological processes and pathways were enriched with $|$normalized enrichment score (NES)$|$ $>$1, false discover rate (FDR) 0.25, and p 0.05, after performing 1000 permutations. In addition, the R package “plyr”, “ggplot2”, “grid” and “gridExtra” in R software was used to summarize the results of the Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathway analysis. The reference file was downloaded from the Molecular Signatures Database (MSigDB) [20].

The univariate and multivariate Cox regression analyses were conducted with Cox proportional hazard regression models to analyze predictive independent factors. Then, a nomogram that included all independent factors predicting overal survival (OS) was developed using R version 4.0.3 (http://www.r-project.org/) with rms packages.

The immune cell infiltration level of each BRCA sample was quantified using CIBERSORT (http://cibersort.stanford.edu/). The relationship between RBPLSig and immune checkpoint inhibitor (ICI) gene expression was determined using the ggstatsplot package and violin plot visualization. To clinically evaluate the model for breast cancer treatment, we calculated the IC${}_{50}$ of therapeutic antitumor drugs, such as lapatinib, gemcitabine, paclitaxel, sunitinib, and tipifarnib, which are recommended for BRCA treatment by AJCC guidelines. The difference in the IC${}_{50}$ between the high- and low-risk groups was compared by a Wilcoxon signed-rank test, and the results are shown as box drawings obtained using pRRophetic and ggplot2 of R. In addition, the consistency between the actual survival and the nomogram was evaluated using calibration curves to predict survival probability. The better the fit between the predicted calibration curve and the standard curve, the better the conformity of the prediction model.

Triple-negative breast cancer (TNBC) cell lines, MDA-MB-231 and HST578T cells, were obtained American Type Culture Collection (ATCC, Manassas, VA, USA) and cultured in DMEM (Gibco, Waltham, MA, USA), containing 4500 mg/L glucose and 10% fetal bovine serum in 5% CO${}_2$ conditions at 37 °C. To verify the function of WAC-AS1, specific siRNA (Guangzhou RiboBio Co., Ltd, Guangzhou, China) was designed and synthesized against it. TRIzol Reagent (Invitrogen, Carlsbad, CA, USA) was used to extract total RNA from two TNBC cell lines; subsequently, the iScript cDNA synthesis kit (Forma Scientific, Marietta, OH, USA) was used to generate cDNA. Real-time PCR was then performed with SYBR Green Real-time PCR Master Mixes (Thermo Fisher Scientific, Waltham, MA, USA) on a CFX96 Real-time PCR Detection System (Bio-Rad, Hercules, CA, USA). Data were normalized to GAPDH as an internal control, and results were analyzed as previously described [21].

The Cell Counting Kit-8 (CCK-8) assay was used to examine the survival of TNBC cells following the protocol of the CCK-8 assay kit (Dojindo, Kyushu, Japan). Survival rates of the two TNBC cell lines after si-WAC-AS1 knockdown was measured at 450 nm [22]. Cellular apoptosis was quantified via terminal deoxynucleotidyl transferase-mediated dUTP nick-end labeling (TUNEL) staining following the protocol of the One Step TUNEL Apoptosis Assay Kit (Keygen, Nanjing, China) [22]. The treated TUNEL-positive TNBC cells were counted under a fluorescence microscope (Olympus, Tokyo, Japan).

We used the cell miner interface (https://discover.nci.nih.gov/cellminer/), containing 727 FDA-approved drugs and clinical trial drugs along with their corresponding protein targets in 60 cancer cell lines, to download the DTP NCI-60 database and RNA-seq data. Then, Spearman or Pearson correlation were used to analyze the association between RBP gene (NSUN6) expression and drug treatment response in five BRCA cell lines.

The 2D structure of the drug molecules was searched using PubChem (https://pubchem.ncbi.nlm.nih.gov/), and a mol2 file was converted using Chem3D (PerkinElmer, Waltham, MA, USA). Receptor proteins were searched from the Protein Data Bank (http://www.rcsb.org/). Dehydration and ligand elimination were performed using PyMOL software (PyMOL2.3, https://pymol.org). The docking site was set in a cubic box in the center of the initial ligand, and a grid map of each atom type in the box was computed. The AutoDockTools 1.5.6 software (https://autodocksuite.scripps.edu/adt/) was used to simulate the molecular docking of potential targets and components. The best scoring conformer of each compound was analyzed and visualized in AutoDockTools-1.5.6 and PYMOL. Concurrently, we used Gromos96 (G43A1, http://gromos.net/) for dynamic simulation, which was performed for 10 ns in combination with the Extended Simple Point Charge (SPCE) water model [23]. System stability was evaluated by the root-mean-square deviation (RMSD) method.

Statistical analyses were conducted using the R software V4.2.1 (http://www.bioconductor.org/) and SPSS V18.0 (IBM Corp., Armonk, NY, USA). The Kaplan–Meier method was used for survival analysis. The statistical analysis of clinical information was performed using a chi-square test or Fisher’s exact test. The statistical significance of the difference between two groups was assessed using a Student’s t-test, and p 0.05 was considered statistically significant.

We first screened 1542 RBP-related encoding genes (mRNAs) and identified the available expression data of 1550 total lncRNAs as RBP-related lncRNAs in the TCGA-BRCA set ($|$R${}^2$$|$ $>$0.3, p 0.001). Among the lncRNAs selected, 61 had prognostic value for patients with BRCA through univariate Cox regression (p 0.05, Supplementary Table 1). Next, 20 prognostic RBP-related lncRNAs were identified as poor prognostic factors though multivariate Cox proportional hazards regression analysis, including MAPT-AS1, AL138724.1, AC004585.1, LINC00667, LINC01871, AL359752.1, AP003469.4, AL122010.1, AC061992.1, AL357054.4, LINC00987, SEMA3B-AS1, WAC-AS1, AL136368.1, AL136531.1, USP30-AS1, Z68871.1, AC245297.3, AC009119.1, and AP003119.3. Supplementary Fig. 1A shows the co-expression network of six lncRNAs (AC004585.1, AP003469.4, WAC-AS1, Z68871.1, AC009119.1, and AP003119.3) and related RBPs. These lncRNAs were identified as poor prognostic factors because HR (hazard ratio) $>$1. The Sankey diagram show the co-expression relationship between RBPs genes and related prognostic lncRNAs Supplementary Fig. 1B).

After multivariate Cox proportional hazards regression analysis, the 20 lncRNAs were used to establish a prognostic RBP-related lncRNA signature (RBPLSig). The estimated RBPLSig was calculated as follows: = (–0.258 $\times$ MAPT-AS1${}_{\mathrm{\text{expression}}}$) + (–0.312 $\times$ AL138724.1${}_{\mathrm{\text{expression}}}$) + (0.248 $\times$ AC004585.1${}_{\mathrm{\text{expression}}}$) + (–0.109 $\times$ LINC00667${}_{\mathrm{\text{expression}}}$) + (–0.265 $\times$ LINC01871${}_{\mathrm{\text{expression}}}$) + (–0.166 $\times$ AL359752.1${}_{\mathrm{\text{expression}}}$) + (0.105 $\times$ AP003469.4${}_{\mathrm{\text{expression}}}$) + (–0.193 $\times$ AL122010.1${}_{\mathrm{\text{expression}}}$) + (–0.231 $\times$ AC061992.1${}_{\mathrm{\text{expression}}}$) + (–0.249 $\times$ AL357054.4${}_{\mathrm{\text{expression}}}$) + (–0.159 $\times$ LINC00987${}_{\mathrm{\text{expression}}}$) + (–0.061 $\times$ SEMA3B-AS1${}_{\mathrm{\text{expression}}}$) + (0.060 $\times$ WAC-AS1${}_{\mathrm{\text{expression}}}$) + (–0.382 $\times$ AL136368.1${}_{\mathrm{\text{expression}}}$) + (–0.320 $\times$ AL136531.1${}_{\mathrm{\text{expression}}}$) + (–0.181 $\times$ USP30-AS1${}_{\mathrm{\text{expression}}}$) + (0.529 $\times$ Z68871.1${}_{\mathrm{\text{expression}}}$) + (–0.109 $\times$ AC245297.3${}_{\mathrm{\text{expression}}}$) + (0.046 $\times$ AC009119.1${}_{\mathrm{\text{expression}}}$) + (0.049 $\times$ AP003119.3${}_{\mathrm{\text{expression}}}$). Using the median risk score as the threshold, the risk score of each patient was calculated, after which patients were divided into low- and high-risk groups (Fig. 1). Notably, patient mortality risk increased with an increasing risk score, and the expression of 20 lncRNAs in the RBPLSig was associated with the risk score (Fig. 1A–C). Kaplan–Meier survival curves showed a significant difference in OS between the high-RBPLSig and low-RBPLSig groups for patients with BRCA (p 0.0001, Fig. 1D), suggesting that the patients with low-risk scores had a greater chance of achieving the same survival time than the high-risk score group (Fig. 1D).

Fig. 1.

The prognostic value of the RBPLSig. (A) Heatmap demonstrating the expression levels of 20 prognostic lncRNAs between the high- and low-risk groups. (B,C) Hierarchical clustering analysis of survival status of patients with increased risk score. (D) Kaplan–Meier analysis of the prognostic model in TCGA-BRCA.

Univariate and multivariate Cox regression analyses were performed to determine whether the RBPLSig is an independent prognostic factor for BRCA. The hazard ratio (HR) of the risk score was 1.104 (95% confidence interval [CI] 1.081–1.128, p 0.001) in the univariate Cox regression analysis (Fig. 2A) and 1.078 (95% CI 1.052–1.104, p 0.001) in the multivariate Cox regression analysis (Fig. 2B). The area under the curve (AUC) values of multiple clinicopathological factors such as age, stage, and Classification of Malignant Tumors (TNM) stage, were significantly different between the two groups stratified by RBPLSig (Fig. 2C). The risk score showed a relatively higher AUC value with respect to multiple clinicopathological factors for the TCGA-BRCA cohort (5-year: 0.777, 8-year: 0.8, 10-year: 0.811, and 15-year: 0.804, Fig. 2D) and validation cohort (GSE58812) (3-year: 0.682, 5-year: 0.719, 8-year: 0.663, and 10-year: 0.657, Fig. 2E). PCA further verified that the RBPLSig was a risk factor related to prognosis, with moderate sensitivity and specificity. The high- and low-risk groups could not be effectively discriminated using the whole genome or RBP-related lncRNAs; however, using RBPLSig, the high- and low-risk patients could be clearly distinguished, further supporting the accuracy of the model. These results indicated that RBPLSig is a significant independent prognostic risk factor for patients with BRCA (Fig. 2F–H).

Fig. 2.

The RBPLSig is an independent prognostic factor for OS prediction. (A,B) Univariate and multivariate Cox regression analyses to screen OS-related factors. (C) Receiver operating characteristic (ROC) curve analysis of clinicopathological features. (D) Time-dependent ROC curves of OS at 5-, 8-,10- and 15-years in TCGA-BRCA cohort. Their AUC values are all greater than 0.75. The closer the AUC is to 1.0, the higher the authenticity of the detection method. (E) Time-dependent ROC curves of OS at 3-, 5-, 8- and 10-years in the validation cohort. Their AUC values are all greater than 0.65. (F–H) Principal component analysis (PCA) of the high- and low-risk groups based on the whole-genome (F), RBP-related lncRNAs (G), and RBPLSig expression profiles (H).

To further assess whether the 20 RBP-related lncRNAs participated in the development of BRCA, the relationship between the risk score and clinicopathological factors was analyzed. As shown in Fig. 3A–E, an association was found between the risk score and clinicopathological factors (p 0.001). Apparently, patients with low-risk scores were clustered more in luminal A and triple negative breast cancer (TNBC) subtypes, whereas patients with high-risk scores were clustered more in HER2-enriched and luminal B subtypes (Fig. 3F, p 0.001).

Fig. 3.

Prognostic risk score associated with the clinical characteristics of patients with BRCA. (A) Heatmap and clinicopathologic features of high- and low-risk groups. *p 0.05, **p 0.01, and ***p 0.001. (B–F) Distribution of risk scores stratified by (B) N-stage (p 0.001), (C) T-stage (p 0.01), (D) stage (p 0.001), (E) age (p 0.05) and (F) tumor subtypes (p 0.001).

GSEA was performed to further examine RBPLSig-associated signaling pathways in BRCA. We observed that transforming growth factor $\beta$ signaling pathways were enriched in the high-risk group. Meanwhile, as shown in Fig. 4A and Supplementary Table 2, several immuno-related signaling pathways, such as B cell receptor, T cell receptor, natural killer (NK) cell-mediated cytotoxicity, and the vascular endothelial growth factor (VEGF) signaling pathway were enriched in the low-risk group. Further, among the lncRNAs selected in RBPLSig, six of these 20 lncRNAs (USP30-AS1, AC004585.1, AC009119.1, AL138724.1, LINC00987 and LINC01871) (Supplementary Fig. 2) were associated with lymph node metastasis. Additionally, 18 of 20 lncRNAs in RBPLSig (AL136368.1, AL136531.1, AL138724.1, LINC01871, MAPT-AS1, AL357054.4, AC061992.1, AL122010.1, USP30-AS1, AL359752.1, LINC00987, AC245297.3, LINC00667, SEMA3B-AS1, AP003119.3, WAC-AS1, AC004585.1 and Z68871.1) were associated with immune subtypes (Fig. 4B, p 0.05). These results indicate that RBPLSig and RBPLSig-associated lncRNAs may participate in tumor progression and the tumor immune microenvironment (TIME). In addition, we analyzed the correlation between 20 lncRNAs and stromal score. The expression level of WAC-AS1, LINC01871, AC004585.1, SEMA3B-AS1, AL136531.1, AP003119.3, AL136368.1, USP30-AS1 and LINC00987 in BRCA patients demonstrated a significant correlation with stromal score (p 0.05), especially AC004585.1 and LINC00987 (r $>$ 0.4) (Supplementary Fig. 3). We also found that eight out of 20 lncRNAs were associated with immune score, and 11 out of 20 lncRNAs were associated with tumor purity (Estimate score) (Supplementary Fig. 3).

Fig. 4.

Bioinformatics analysis between the two groups stratified by RBPLSig. (A) GSEA for biological pathways and processes correlated with immune score values in the cohort from TCGA. (B) Association of 20 lncRNAs expression with different immune infiltrate subtypes tested with ANOVA (analysis of variance) in TCGA-BRCA. ns, not significant; *p 0.05; **p 0.01; ***p 0.001. C1: wound healing, C2: INF-r dominant, C3: inflammatory, C4: lymphocyte depleted, C5: immunologically quiet, and C6: TGF-$\beta$-dominant.

To evaluate the correlation between the risk score and tumor immune cell infiltration, the CIBERSORT algorithm was used to investigate the scale of 22 immune cell types. Patients with high-risk scores had higher macrophage M0 and macrophage M2 infiltration levels than those with low-risk scores (p $\le$ 0.001). In contrast, patients with low-risk scores had higher ratios of B cell-naive, plasma cells, T cells CD8, T cells CD4 memory activated, T cells regulatory, NK cell activated, dendritic cells resting, and neutrophils (Fig. 5A, p $\le$ 0.05). This indicates that patients with different immune status patterns in different groups could be classified by the RBPLSig.

Fig. 5.

Different responses to immunotherapy, chemotherapy and targeted therapy between the two groups stratified by RBPLSig. (A) Violin plot of 22 immune-related terms were incorporated to assess the abundance of immune cells between the two groups stratified by RBPLSig, in which red represents the high-risk samples, and blue the low-risk samples. (B) The expression of PD-1 (p 0.001), PD-L1 (p 0.001), CTLA-4 (p 0.001) and Siglec-15 (p 0.05) between the low- and high-risk groups in the cohort from TCGA. (C) Five drugs include lapatinib (p 0.001), gemcitabine (p 0.001), paclitaxel (p 0.001), sunitinib (p 0.001) and tipifarnib (p 0.001) susceptibility distribution between the high and low risk groups.

Expression of the ICI (immune checkpoint inhibitor) genes and/or their ligands have attracted widespread attention as BRCA biomarkers for immune checkpoint blockade therapy. It is critically important to effectively screen BRCA patients who might benefit most from different ICI immunotherapies. To further investigate the association between risk score and ICIs, cytotoxic T-lymphocyte-associated antigen 4 (CTLA-4), sialic acid binding Ig-like lectin 15 (Siglec-15), programmed cell death-1 (PD-1), and programmed cell death-ligand 1 (PD-L1) were selected to compare the expression patterns between the two groups stratified by the RBPLSig. The result indicates that patients in the low-risk group tended to have higher PD1, PD-L1, and CTLA-4 expression levels than those in the high-risk group, whereas patients in the high-risk group tended to have higher Siglec-15 expression levels (Fig. 5B, p 0.05). These results further indicate that the RBPLSig may serve as a potential predictive biomarker for treatment responses to ICI immunotherapy. We validated the RBPLSig for 12 cancer types suitable for immunotherapy to determine whether it could be effectively used as a prognostic signature for other cancer types. Interestingly, of the 12 cancer types, RBPLSig was negatively correlated with the prognosis of 11 cancer types, while the AUCs corresponding to the 1-, 3-, and 5-year survival for five cancers (CESC: cervical squamous cell carcinoma and endocervical adenocarcinoma; COAD: colon adenocarcinoma; KIRC: kidney renal clear cell carcinoma; KIRP: kidney renal papillary cell carcinoma and LAML: acute myeloid leukemia) were higher than 0.7 (Table 1). In addition to immune checkpoint blockade therapy, the associations between the risk score and the efficacy of targeted therapeutics and chemotherapeutics in treating patients with BRCA were then identified. The low-risk group had a higher IC${}_{50}$ for lapatinib (p 0.001). However, the high-risk group had a higher IC${}_{50}$ for gemcitabine, paclitaxel, sunitinib, and tipifarnib (p 0.001). Collectively, this suggests that the RBPLSig may be a sensitive predictor of drug sensitivity (Fig. 5C).

To generate a clinically applicable model for individual OS prediction, we first constructed a nomogram that accurately estimates the 5-, 8-, 10-, and 15-year survival probabilities in BRCA patients using the RBPLSig and other clinicopathological independent prognostic factors (Fig. 6A). The calibration curves of the nomogram for predicting patient survival at 5-,8-,10-, and 15-years are shown in Fig. 6B–E. Considering that (1) the risk score was lower for luminal A and TNBC than for HER2-enriched and luminal B subtypes, and (2) patients with a low-risk score had an improved TIME (appropriate for immunotherapy), we further constructed a nomogram for BRCA subtypes (luminal A and TNBC) appropriate for immunotherapy evaluation (Supplementary Fig. 4). To further explore the value of predictive nomograms, we developed an online analysis tools to predict the survival rate and recommended medication for BRCA patients (https://www.origingenetic.com/BreastCancerModel). With the aid of the prediction online tools, clinicians can assist in evaluating patient survival rates to optimize treatment plans.

Fig. 6.

Construction of a nomogram for survival prediction. (A) The clinical prognostic nomogram developed to predict 5-, 8-, 10-, and 15-year survival. (B–E) Calibration curves showing nomogram predictions for (B) 5-, (C) 8-, (D) 10- and (E) 15-year survival.

For further study on the potential function of these lncRNAs and to verify the aforementioned results, we focused on the upregulated WAC-AS1, which is associated with metastasis and poor prognosis for OS in selective patient subgroups from TCGA and Gene Expression Omnibus (GEO) datasets (Fig. 7A–C).

Fig. 7.

WAC-AS1 was upregulated in TNBC tissues and WAC-AS1 silencing induced TNBC cells apoptosis. (A) Microarray analysis of WAC-AS1 in TCGA-TNBC tissues corresponding normal tissues. (B,C) Kaplan-Meier survival curves of patients with different level of WAC-AS1 in TCGA (B) and GEO (C) datasets. (D) TNBC cell lines (MDA-MB-231and HS578T) was transfected with si-WAC-AS1. The expression of WAC-AS1 in TNBC cell lines was analyzed by qRT-PCR. (E–G) The Effects of si-WAC-AS1 on the proliferation and apoptosis of TNBC cell lines was detected by CCK-8 assay (E,F) and TUNEL assay (G). * p 0.05, **p 0.01 compared with si-NT group. (H,I) GSEA was used to analyze the signaling pathways enrichment in different groups.

Next, we synthesized a specific smart RNA silencer (including small interfering RNA and small nucleolar RNAs) against WAC-AS1. The RT-qPCR indicates that smart RNA silencer-WAC-AS1 (si-WAC-AS1) could abolish the expression of WAC-AS1 in TNBC cell lines (MDA-MB-231 and HST578T cells; Fig. 7D). The CCK-8 assay indicates that cell viability was diminished by silencer-WAC-AS1 transfection compared to that with si-NT (Fig. 7E,F, p 0.05). The TUNEL assay was in line with the CCK-8 assay data and demonstrated that MDA-MB-231 and HST578T cell apoptosis could be enhanced by si-WAC-AS1 transfection compared to that with si-NT (Fig. 7G, p 0.05). Finally, in the pathway analysis, different malignancy-associated pathways were enriched in different groups based on WAC-AS1 (Fig. 7H,I and Supplementary Table 3).

Currently, lncRNAs are usually used as prognostic or diagnostic biomarkers and lncRNA-based drug targets [24], the latter of which, while reported in several studies, still requires more research [25]. Accordingly, the search for new lncRNA-based drugs is considered a promising field. Targeting the binding partners (pro-oncogenic RBPs) may indirectly affect the function of these related pro-oncogenic lncRNAs. Drug sensitivity analysis was performed on five breast cancer cell lines, including MCF7, MDA-MB-231, HS 578T, BT-549, and T-47D ($|$r$|$ $\ge$0.9, p 0.001). Strikingly, after further screening, NOP2/Sun RNA methyltransferase 6 (NSUN6; PDBID: 5WWQ), which was co-expressed with the pro-oncogenic lncRNA WAC-AS1 (cor = 0.37, p 0.001), was associated with sensitivity to navitoclax in five BRCA cell lines (cor = 0.996, p 0.001, Supplementary Fig. 5). We gained mechanical insights into the potential mechanisms of the druggable compounds of an RBP gene, putative NSUN6, using molecular docking. The crystal structure of the RBP gene, mainly protease, was collected from the PDB database (http://www.rcsb.org/structure/NSUN6) with PDB ID 5WWQ, for docking analysis with navitoclax (Fig. 8). The active cavity box parameter setting center_x, y, z was 54.315, 33.779 and 135.063, respectively, while size_x, y, z was 92, 64, 108, whereas amino acid residues were within 3.4Å around the active site. The minimum docking affinity between navitoclax and NSUN6 is –9.0 (kcal/mol). The lower the docking score (the greater the negative value), the higher the binding force is between the compound and the protein. The bond action showed that navitoclax forms a hydrogen bond with the residue of six amino acids (SER-210, SER-207, TYR-131, ARG-197, GLN-78, and GLU-74) of NSUN6(RBP) (Fig. 8A), suggesting that navitoclax has a good binding activity with the NSUN6 protein. Next, the behavior of NSUN6, associated with navitoclax sensitivity in a simulative aquatic ecosystem, was analyzed. After a 10 ns simulation, the root-mean-square deviation (RMSD) tended to be flat; the display system reached equilibrium. Although the structure and conformation of the ligand changed, the complex did not separate, and the ligand receptor was tightly bound in the target region, indicating that the site changed to be a possible target (Fig. 8B,C). Taken together, these predictions indicate that navitoclax may be a useful treatment option for TNBC patients with high expression of the pro-oncogenic NSUN6 and the binding partner lncRNAs. Although preliminary, this study presents compelling data for the design of lncRNA-based drugs; it is hoped that future studies will expand upon our themes.

Fig. 8.

Molecular docking analysis. (A) Molecular docking results of navitoclax-NSUN6 (NOP2/Sun RNA methyltransferase 6 NSUN6; PDBID: 5WWQ). Hydrogen bonds are shown as yellow dotted lines. (B) Snapshots of navitoclax-NSUN6 system after 10 ns. (C) RMSD of navitoclax-NSUN6 systems. RBP, RNA-binding protein.