A Risk Score System Based on the Methylation Levels of 15 RNAs in Breast Cancer

Abstract

Background:

Breast cancer (BC) occurs in the epithelial tissues of the breast gland, which is the most common cancer in women. This study is implemented to construct a risk score system for BC.

Materials and Methods:

The methylation data of BC from The Cancer Genome Atlas database (the training set) and GSE37754 from Gene Expression Omnibus database (the validation set) were downloaded. The differentially methylated RNAs (DMRs) between BC and normal samples were screened by limma package, and the correlations between the expression levels and methylation levels of the DMRs were analyzed to calculate their Pearson correlation coefficients (PCCs) using the cor.test function. To build the risk score system, the optimal RNAs were identified by penalized package. Subsequently, the nomogram survival model was established using the rms package. The lncRNA-mRNA comethylation network was constructed by Cytoscape software, and then enrichment analysis was performed using DAVID tool.

Results:

From the 1170 DMRs between BC and normal samples, 800 DMRs with significant negative PCCs were screened. For building the risk score system, the 15 optimal RNAs were selected. Afterward, the nomogram survival model based on four independent clinical prognostic factors (including age, radiation therapy, tumor recurrence, and RS model status) was constructed. In the comethylation network, the long noncoding RNA (lncRNA) PRNT was comethylated with FAM19A5 and RBM24. For the mRNAs in the comethylation network, angiogenesis and pathways in cancer were enriched.

Conclusion:

The risk score system and the nomogram survival model might be of great importance for the prognosis prediction of BC patients.

Background

Breast cancer (BC) is the most common cancer in women, which has the symptoms of breast lump, breast deformation, skin dimpling, and nipple discharge.¹ The main risk factors of BC include family history, ionization radiation, obesity, lack of exercise, excessive alcohol and cigarette consumption, early menstruation, and later birth or infertility.² BC takes up 25% of all tumor cases in women, which was responsible for ∼459,000 deaths worldwide in 2015.^3,4 The curative effects of BC are mainly related to tumor stage and onset age, and early diagnosis is an important way to improve prognosis and reduce mortality of BC.⁵ Therefore, more studies should be conducted to promote the diagnosis and prognosis of BC.

DNA methylation can promote carcinogenesis and tumor progression through inducing abnormal gene expression and genomic instability, and methylation detection is valuable for screening tumor susceptible people.^6,7 The promoter methylation of cadherin 1 (CDH1) is related to estrogen receptor status and lymph node metastasis of BC patients, which independently predicts unfavorable outcome in the tumor.⁸ Ras association domain family member 2 (RASSF2) hypermethylation in luminal tumors is more frequent and is correlated with favorable prognosis, and thus methylation status of RASSF2 can serve as a marker for prognosis in BC.^9,10 The promoter methylation of the long noncoding RNA (lncRNA) long intergenic nonprotein coding RNA 472 (LINC00472) has an influence on LINC00472 expression in BC, and its expression indicates the recurrence of grade 2 BC patients.¹¹ The lncRNA maternally expressed 3 (MEG3) acts as a tumor suppressor in the metastasis of BC, and the miR-506/SP3/SP1/DNA methyltransferase 1 (DNMT1)/MEG3 axis plays roles in mediating the invasion and migration of BC cells.¹² Nevertheless, the roles of lncRNA and mRNA methylation in the prognosis of BC have not been fully understood.

The polygenic risk score involving multiple single nucleotide polymorphisms can improve the predictive ability of the Breast Cancer Surveillance Consortium risk model, which contributes to select the appropriate treatment modalities for BC patients.^13
–15 However, the risk score system based on the methylation levels of lncRNAs and mRNAs has not been established. In this study, the differentially methylated RNAs (DMRs) between BC and normal samples in The Cancer Genome Atlas (TCGA) dataset were explored. Then, the independent prognosis-associated DMRs were further screened for constructing the risk score system. Moreover, nomogram survival model was built based on the independent clinical prognostic factors. In addition, lncRNA-mRNA comethylation network analysis and enrichment analysis were conducted to explore the potential functions of the RNAs. The findings from this study might help to improve the accuracy of prognostic prediction in BC patients.

Materials and Methods

Data source

From TCGA (https://cancergenome.nih.gov/) database, the methylation data of BC (platform: Illumina Infinium Human Methylation 450 BeadChip; downloaded in September 10, 2019; including 890 samples) was extracted. After corresponding to the clinical information provided for the samples, 872 samples (including 776 BC samples with survival prognosis information and 96 normal samples) were included into this study as a training set. Meanwhile, the methylation data under GSE37754 (platform: Illumina Infinium Human Methylation 450 BeadChip) was obtained from Gene Expression Omnibus (GEO, www.ncbi.nlm.nih.gov/geo/) database. Among GSE37754, there were 62 BC samples with survival prognosis information and 10 normal samples (the validation set).

Differential analysis

The methylated sites in the TCGA dataset were annotated according to the platform annotation information, and only the methylated sites involving specific genes and lncRNAs were retained. Using the R package limma (version 3.34.7),¹⁶ the DMRs were screened from the BC samples and normal samples. The false discovery rate (FDR) <0.05 and |log₂ fold change (FC)| > 1 were selected as the thresholds.

Correlation analysis

The BC samples with both methylation data and expression data were selected, and then the expression levels and methylation levels of the identified DMRs were obtained. Using the cor.test function¹⁷ in R, correlation analysis was conducted for the overall expression levels and methylation levels of the DMRs to calculate their Pearson correlation coefficients (PCCs). After that, PCC was calculated for the expression level and methylation level of a single DMR. Finally, only the DMRs with significant negative PCCs were retained for further analysis.

Construction of risk score system

Combined with the clinical prognosis information of the 776 BC samples in the training set, the DMRs with methylation levels significantly correlated with prognosis were screened using the univariate Cox regression analysis in the R package survival.¹⁸ Afterward, the independent prognosis-associated DMRs were further selected using the multivariate Cox regression analysis in the survival package.¹⁸ The log-rank p-value <0.05 was selected as the threshold of significant correlation.

Using the Cox-Proportional Hazards (Cox-PH) model in the R package penalized¹⁹ (version 0.9.50), the optimal combination of the independent prognosis-associated DMRs was selected based on their methylation levels. The optimized parameter “lambda” in the model was obtained through 1000 times of cycle calculation using crossvalidation likelihood (cvl).²⁰

On the basis of the prognostic coefficients and methylation levels of the DMRs in the optimal combination, the following risk score system was established:

Here, β_DMRs and Methylation _DMRs separately represent the prognostic coefficients and methylation levels of the optimal DMRs in the training set.

The RS of each sample in the training set was calculated, and then the samples were classified into high-risk and low-risk groups based on the RS median value. Using the Kaplan–Meier (KM) method in survival package,¹⁸ the correlation between high- and low-risk grouping and actual survival prognosis information was evaluated. Similarly, the methylation levels of the optimal DMRs were extracted from the samples in the validation set, and the samples were divided into high- and low-risk groups according to their RS median value. Finally, the correlation between the high- and low-risk grouping and the actual survival prognosis information of the samples in the validation set was assessed using KM method in survival package.¹⁸

Construction of nomogram survival model

Using the univariate and multivariate Cox regression analyses in survival package,¹⁸ the independent clinical prognostic factors in the training set were screened. The threshold was defined as the log-rank p-value <0.05. To explore the correlation between the independent clinical prognostic factors and survival prognosis, the R package rms (version 5.1–2)²¹ was used to build the nomogram survival model combined with the risk information discriminated by the risk score system.

Moreover, the model based on the independent clinical prognostic factors and the model based on both the independent clinical prognostic factors and the optimal DMRs were constructed. In addition, area under the receiver operating characteristic (ROC) curve (AUROC) based on pROC package²² (version 1.14.0) and C-index based on survcomp package²³ (version 1.34.0) were calculated to compare the efficiencies of the two kinds of models with that of the risk score system.

LncRNA-mRNA comethylation network analysis and enrichment analysis

The methylation levels of the DMRs significantly correlated with prognosis were extracted from the training set, and then the PCCs of the methylation levels of the lncRNAs and mRNAs among these DMRs were calculated by the cor.test function.¹⁷ The lncRNA-mRNA pairs with |PCC| > 0.4 were selected. Using Cytoscape software²⁴ (version 3.6.1), lncRNA-mRNA comethylation network was constructed. Under the cutoff criterion of p-value <0.05, Gene Ontology (GO)_Biology Process (BP) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses for the mRNAs in the comethylation network were performed using DAVID tool²⁵ (version 6.8, https://david.ncifcrf.gov/).

Results

Differential analysis

Based on the platform annotation information, a total of 14,005 methylated sites involving specific genes and lncRNAs were obtained. There were 1170 DMGs (including 16 lncRNAs and 1154 mRNAs) between BC samples and normal samples (Fig. 1A, B). After the top 50 hypermethylated and hypomethylated RNAs were selected according to the ascending order of log₂FC, the bidirectional hierarchical clustering heatmap was drawn based on their methylation levels. The clustering heatmap indicated that the samples were clearly divided into two different groups (Fig. 1C).

FIG. 1.

The results of differential analysis. (A) The Volcano plot for the DMRs (Gray and black dots separately represent DMRs and non-DMRs; the horizontal dashed line indicates FDR <0.05, and the vertical dashed lines represent |log₂ fold change (FC)| > 1); (B) The distribution of different kinds of DMRs, including lncRNAs and mRNAs; (C) The clustering heatmap based on the methylation levels of the top 50 hypermethylated and hypomethylated RNAs. Black and white sample strips separately represent tumor samples and normal samples. The different color represents different methylated level. DMR, differentially methylated RNA; FDR, false discovery rate; lncRNA, long noncoding RNA.

Correlation analysis

Correlation analysis showed that the overall expression levels and methylation levels of the DMRs had significant negative correlation (Cor. = −0.2549, p-value <2.20e-16) (Fig. 2). Through calculating the PCC for the expression level and methylation level of a single DMR, 800 DMRs with significant negative PCCs were selected for the subsequent analysis.

FIG. 2.

The scatter diagram showing the correlation between the overall expression levels and methylation levels of the DMRs in the training set. The gray dashed indicates the distribution trend of the dots.

Construction of risk score system

Combined with univariate Cox regression analysis, 175 DMRs significantly correlated with prognosis were selected. From the 175 DMRs, 31 independent prognosis-associated DMRs (including four lncRNAs and 27 mRNAs) were further screened by multivariate Cox regression analysis.

Based on the 31 DMRs, the optimal combination involving 15 independent prognosis-associated DMRs (one lncRNAs and 14 mRNAs) was identified using the Cox-PH model (protective RNAs: prion locus lncRNA, testis expressed, PRNT; tuftelin 1, TUFT1; zinc finger protein 462, ZNF462; chromobox 5, CBX5; solute carrier family 25 member 33, SLC25A33; and TOX; risky RNAs: acetyl-CoA carboxylase beta, ACACB; glutathione S-transferase mu 1, GSTM1; family with sequence similarity 19 (chemokine (C-C motif)-like), member A5, FAM19A5; zinc finger protein 804B, ZNF804B; RNA-binding motif protein 24, RBM24; potassium voltage-gated channel modifier subfamily V member 1, KCNV1; sarcoglycan beta, SGCB; cysteine-rich secretory protein LCCL domain-containing 1, CRISPLD1; and cell adhesion molecule L1 like, CHL1) (Fig. 3A and Table 1). Then, the risk score system based on the prognostic coefficients and methylation levels of the optimal DMRs was built. According to the RS median values, the samples in the training set and validation set were classified into high-risk group and low-risk group. The RSs, overall survival time, and the methylation levels of the optimal DMRs are also presented (Fig. 3B, C).

FIG. 3.

The coefficients and methylation levels of the optimal RNAs, and the risk scores and OS time of the samples in the training set. (A) The coefficient curve for the optimal RNAs; (B) The distribution diagrams of the risk scores (top) and the OS time (middle) of the samples in TCGA, and the heatmap for the methylation levels of the optimal RNAs (bottom). (C) The distribution diagrams of the risk scores (top) and the OS time (middle) of the samples in GSE37754. TCGA, The Cancer Genome Atlas; OS, overall survival.

Table 1.

The Differentially Methylated RNAs Involved in the Optimal Combination

Symbol	ID	Chromosome location		Type	Univariate Cox regression analysis			LASSO coefficient
Symbol	ID	Chr.	Position	Type	HR	95% CI	p-value	LASSO coefficient
Protective RNAs
PRNT	cg26097573	chr20	4669490	lncRNA	0.678	0.257–0.788	4.32E-02	−0.7041846
TUFT1	cg27577972	chr1	149780126	Coding	0.0218	0.00139–0.341	6.43E-03	−1.2737752
ZNF462	cg13804575	chr9	108666820	Coding	0.315	0.148–0.672	2.53E-03	−1.2534855
CBX5	cg27408285	chr12	52939631	Coding	0.368	0.168–0.806	1.18E-02	−1.0428817
SLC25A33	cg26607031	chr1	9525288	Coding	0.241	0.0938–0.623	3.09E-03	−0.6746042
TOX	cg27496650	chr8	60194800	Coding	2.459	1.103–5.480	2.69E-02	−0.5386747
Risky RNAs
ACACB	cg24284369	chr12	108077370	Coding	4.907	1.467–16.41	8.23E-03	0.4775332
GSTM1	cg24506221	chr1	110031924	Coding	1.836	1.003–3.362	4.78E-02	0.6070588
FAM19A5	cg27652464	chr22	47265439	Coding	5.403	1.924–15.18	1.20E-03	0.6176417
ZNF804B	cg27319898	chr7	88226939	Coding	13.09	4.597–37.25	6.28E-07	0.6212348
RBM24	cg25302957	chr6	17390263	Coding	3.644	1.309–10.14	1.21E-02	0.6958808
KCNV1	cg23291534	chr8	111057219	Coding	3.315	1.632–6.734	7.86E-04	1.0075978
SGCB	cg26148714	chr4	52599332	Coding	21.15	5.887–75.99	3.08E-07	1.0321157
CRISPLD1	cg21844856	chr8	76059200	Coding	31.26	6.367–153.4	1.27E-05	1.7010498
CHL1	cg25482786	chr3	213496	Coding	9.409	3.827–23.13	5.09E-08	1.8022721

ID, identification; HR, Hazard Ratio; CI, confidence interval.

Combined with KM method, the correlation between the risk grouping and actual survival prognosis information was evaluated in both the training set and the validation set. KM curves suggested that the risk groups divided by the risk score system and actual survival prognosis information had significant correlation in both the training set (Fig. 4A) and the validation set (Fig. 4B).

FIG. 4.

The KM curves and ROC curves. (A) The KM curves (black and gray curves represent low- and high-risk groups, respectively) and ROC curve (black, red and blue separately represent 1-year, 3-year, and 5-year) for the training set; (B) The KM curves (black and red curves represent low- and high-risk groups, respectively) and ROC curve for the validation set. HR, hazard ratio; AUC, area under the receiver operating characteristic curve; KM, Kaplan–Meier; ROC, receiver operating characteristic.

Construction of nomogram survival model

Based on Cox regression analysis, four independent clinical prognostic factors (including age, radiation therapy, tumor recurrence, and RS model status) in the training set were identified (Table 2). The KM curves for age, tumor recurrence, and radiation therapy are shown in Figure 5. The prognosis of the BC patients with lower age, without tumor recurrence, and receiving radiation therapy was better, which was consistent with the actual situation.

FIG. 5.

The KM curves. (A) The KM curves showing the correlation of age with prognosis (black and gray separately represent the patients aged no more than 60 years and over 60 years); (B) The KM curves showing the correlation of recurrence with prognosis (black and gray separately represent without recurrence and with recurrence); (C) The KM curves showing the correlation of radiotherapy with prognosis (black and gray represent without radiotherapy and with radiotherapy, respectively).

Table 2.

The Screening Results of Independent Clinical Prognostic Factors

Clinical characteristics	TCGA (N = 776)	Univariables Cox			Multivariables Cox
Clinical characteristics	TCGA (N = 776)	HR	95%CI	p	HR	95%CI	p
Age(years, mean ± SD)	58.06 ± 13.16	1.033	1.018–1.048	1.48E-05	1.038	1.015–1.061	1.27E-03
ER status(Positive/Negative/−)	566/169/41	0.655	0.425–1.009	6.15E-02	—	—	—
PR status(Positive/Negative/−)	493/239/44	0.762	0.504–1.154	1.99E-01	—	—	—
Her2 status(Positive/Negative/−)	58/404/314	0.929	0.422–2.043	8.54E-01	—	—	—
Histological type(Basal/Her2/LumA/LumB/Normal/−)	87/32/277/126/18/236	1.171	0.939–1.461	1.61E-01	—	—	—
Pathologic_M(M0/M1/−)	609/13/154	5.354	2.715–10.56	8.58E-08	1.506	0.478–4.742	4.84E-01
Pathologic_N(N0/N1/N2/N3/−)	348/265/95/56/12	1.598	1.293–1.974	1.01E-05	1.071	0.708–1.619	7.46E-01
Pathologic_T(T1/T2/T3/T4/−)	199/443/108/23/3	1.504	1.183–1.912	8.03E-04	1.179	0.759–1.833	4.63E-01
Pathologic_stage(I/II/III/IV/−)	156/434/196/11/9	2.321	1.721–3.126	2.21E-08	2.523	1.091–5.830	5.04E-02
Radiotherapy(Yes/No/−)	408/312/56	0.567	0.358–0.899	1.45E-02	0.422	0.234–0.760	4.10E-03
Recurrence(Yes/No/−)	74/616/86	9.109	5.531–15.00	2.00E-16	8.209	4.604–14.64	9.75E-13
Risk status(High/Low)	388/388	3.677	2.339–5.779	8.92E-10	5.484	2.762–10.89	1.16E-06
Dead(Death/Alive)	102/674	—	—	—	—	—	—
Overall survival time(months, mean ± SD)	42.02 ± 38.74	—	—	—	—	—	—

TCGA, The Cancer Genome Atlas; HR, hazard ratio; CI, confidence interval.

Besides, a nomogram survival model was constructed based on the independent clinical prognostic factors (Fig. 6). In the nomogram survival model, the “Total points” axis was used to integrate various clinical indicators to perform survival prediction for the samples. Each independent clinical prognostic factor could be aligned to a point and the total points equal the sum of points of each clinical prognostic factor. The 3-year survival probability and 5-year probability could be predicted based on the total points of each prognostic factor. For example, the total points of a woman 72 years of age (point = 27), did not receive radiation (point = 0), suffered from recurrence (point = 0) and with a RS of 1 (point = 33) is 61. Her 3-year survival probability is around 0.54 and her 5-year survival probability is around 0.37. Furthermore, the predicted 3-year and 5-year survival rates based on the nomogram survival model was compared with the actual survival rates to obtain their consistencies (Fig. 7A). In addition, the efficiencies of the model based on the independent clinical prognostic factors and the model based on both the independent clinical prognostic factors and the optimal DMRs were compared with that of the risk score system. The results showed that the efficiency of the model based on both the independent clinical prognostic factors and the optimal DMRs was the best (AUROC = 0.966, C-index = 0.881) (Fig. 7B).

FIG. 6.

The nomogram survival model based on the independent clinical prognostic factors, including age, radiation therapy, tumor recurrence, and RS model status. Each independent clinical prognostic factor could be aligned to a point and the total points equal the sum of points of each clinical prognostic factor. The 3-year survival probability and 5-year survival probability could be predicted based on the total points. The vertical line illustrates the correspondence between the total points and the survival probability. RS, risk score.

FIG. 7.

The consistency diagrams of survival rates and the comparison of different kinds of models. (A) The consistency diagrams for the predicted 3-year (left) and 5-year (right) survival rates compared with the actual survival rates (the horizontal and vertical axes stand for the predicted OS and the actual OS, respectively); (B) The ROC curves (top) and the comparison parameters (bottom) of different kinds of models. AUROC, area under the receiver operating characteristic curve.

LncRNA-mRNA comethylation network analysis and enrichment analysis

After 157 eligible lncRNA-mRNA pairs were obtained, the lncRNA-mRNA comethylation network (involving 98 nodes and 157 edges) was constructed (Fig. 8). In the network, the lncRNA PRNT was comethylated with FAM19A5 and RBM24. For the mRNAs in the comethylation network, 15 GO_BP terms (such as positive regulation of transcription from RNA polymerase II promoter, p-value = 5.84E-06; adult locomotory behavior, p-value = 6.00E-06; and angiogenesis, p-value = 5.21E-03) and three KEGG pathways (such as neuroactive ligand–receptor interaction, p-value = 3.36E-04; calcium signaling pathway, p-value = 1.12E-02; and pathways in cancer, p-value = 4.21E-02) were enriched (Table 3).

FIG. 8.

The lncRNA-mRNA comethylation network. Squares and circles represent lncRNAs and mRNAs, respectively. Fold change represents the changes of methylation level. The enlarged nodes represent the optimal RNAs.

Table 3.

The Biological Process Terms and Pathways Enriched for the mRNAs in the Comethylation Network

Category	Term	Count	p-value
Biological process	GO:0045944∼positive regulation of transcription from RNA polymerase II promoter	18	5.84E-06
	GO:0008344∼adult locomotory behavior	6	6.00E-06
	GO:0045893∼positive regulation of transcription, DNA-templated	10	1.07E-03
	GO:0019221∼cytokine-mediated signaling pathway	5	4.22E-03
	GO:0043524∼negative regulation of neuron apoptotic process	5	4.33E-03
	GO:0007204∼positive regulation of cytosolic calcium ion concentration	5	4.57E-03
	GO:0001525∼angiogenesis	6	5.21E-03
	GO:0030154∼cell differentiation	8	8.37E-03
	GO:0008284∼positive regulation of cell proliferation	8	8.75E-03
	GO:0045892∼negative regulation of transcription, DNA-templated	8	1.24E-02
	GO:0042493∼response to drug	6	1.82E-02
	GO:0000122∼negative regulation of transcription from RNA polymerase II promoter	9	2.74E-02
	GO:0007268∼chemical synaptic transmission	5	3.23E-02
	GO:0007186∼ G-protein-coupled receptor signaling pathway	10	3.54E-02
	GO:0007267∼ cell–cell signaling	5	3.85E-02
Pathway	hsa04080:Neuroactive ligand–receptor interaction	8	3.36E-04
	hsa04020:Calcium signaling pathway	5	1.12E-02
	hsa05200:Pathways in cancer	6	4.21E-02

GO, gene ontology.

Discussion

In this study, 1170 DMGs between BC samples and normal samples were screened, from which 800 DMRs with significant negative PCCs were selected. After the optimal combination involving 15 independent prognosis-associated DMRs (protective RNAs: PRNT, TUFT1, ZNF462, CBX5, SLC25A33, and TOX; risky RNAs: ACACB, GSTM1, FAM19A5, ZNF804B, RBM24, KCNV1, SGCB, CRISPLD1, and CHL1) was identified, the risk score system was built. Based on the four independent clinical prognostic factors (including age, radiation therapy, tumor recurrence, and RS model status), the nomogram survival model was established. Moreover, the efficiency of the model based on both the independent clinical prognostic factors and the optimal DMRs was better. For the mRNAs in the comethylation network, 15 GO_BP terms (such as angiogenesis) and three KEGG pathways (such as pathways in cancer) were enriched. Uncontrolled angiogenesis is implicated in tumor growth, invasion, and metastasis, and mast cells may affect the adjuvant therapy of BC through inhibiting angiogenesis.²⁶ Thus, angiogenesis and pathways in cancer might be implicated in the prognosis of BC patients.

Although FAM19A5, RBM24, ZNF462, and KCNV1 have not been reported to be correlated with BC, they play roles in other tumors. FAM19A5 is differentially expressed between cholangiocarcinoma (CCA) and benign biliary tract diseases (BBTDs), which may be novel serum marker for distinguishing CCA from BBTDs.²⁷ Through increasing miR-25 expression, RBM24 functions as a tumor suppressor gene in nasopharyngeal carcinoma.²⁸ The lncRNA PRNT was comethylated with FAM19A5 and RBM24 in the lncRNA-mRNA comethylation network, indicating that PRNT might act in BC through comethylating with FAM19A5 and RBM24. ZNF462 is a potential target of miR-210 in pancreatic cancer (PC), and miR-210 expression can be promoted by hypoxia through the HIF-1 alpha-dependent pathway in PC cells.²⁹ KCNV1 is a diagnostic signature for the early screening of bladder cancer, and the methylation of LIM homeobox transcription factor 1 alpha (LMX1A) and ventral anterior homeobox 1 (VAX1) is correlated with the recurrence of the tumor.³⁰ Therefore, the authors speculated that FAM19A5, RBM24, ZNF462, and KCNV1 might affect the prognosis of BC.

The coexpression of TUFT1 and Rab5 indicates poor outcome in triple-negative breast cancer (TNBC) patients receiving chemotherapy, and the TUFT1/Rab5/Rac1 axis may provide promising targets for the therapy of TNBC.^31,32 CBX5 transcriptional and protein levels are reduced in metastatic BC compared with nonmetastatic BC, and the negative relationship between STET and HP1α coding CBX5 mRNA levels is detected in BC.^33,34 Solute carrier family 25A member 43 (SLC25A43) is found to be a putative tumor suppressor gene in epidermal growth factor receptor 2 (HER2)-positive BC, and its methylation is involved in gene silencing in the absence of loss of heterozygosity.³⁵ TOX high-mobility group box family member 3 (TOX3) rs3803662 polymorphism has a significant correlation with the BC risk, which can increase the risk of the tumor.³⁶ These suggested that the protective RNAs, TUFT1, CBX5, SLC25A33, and TOX, might have influences on the prognosis of BC.

ACACB, cyclin B2 (CCNB2), cyclin-dependent kinase 1 (CDK1), mitotic arrest-deficient 2 like 1 (MAD2L1), and peroxisome proliferator-activated receptor gamma (PPARG) are possible diagnosis markers in breast ductal carcinoma, and the complementarity between CDK1 and aminophenazone, rosoxacin, and pomalidomide induces novel pharmacological effects in the tumor.³⁷ Glutathione S-transferase theta 1 (GSTT1)/GSTM1 double null and GSTM1 null polymorphisms in BC may be related to an enhanced tumor response, whereas the GSTM1 null genotype may be correlated with a decreased overall survival.³⁸ CHL1 hypermethylation is detected in BC, which is correlated with downregulated CHL1 expression and shorter progression-free survival in the patients with BC.³⁹ Thus, the risky RNAs ACACB, GSTM1, and CHL1 might be correlated with the prognosis of BC.

Conclusion

In summary, 1170 DMGs between BC samples and normal samples were identified. Besides, the 15 RNA-based risk score system and the nomogram survival model might be valuable for predicting the prognosis of BC patients. Nevertheless, further experiments should be designed and implemented to support these results.

Ethics Approval and Consent to Participate

This study was approved by the Ethics Committee of Beijing Shijitan Hospital, Capital Medical University.

Availability of Data and Materials

Not applicable. This study was only the primary research, and further study has been in progress.

Footnotes

Authors' Contributions

Conception of the research and drafting the article: Y.S. and R.W.; Acquisition of data: Y.S and R.W.; Analysis and interpretation of data: Y.S. and R.W.; Statistical analysis: Y.S. and R.W.; Revision of article for important intellectual content: Y.S. and R.W.

Disclosure Statement

The authors declare that they have no competing interests.

Funding Information

No funding was received for this article.

References

Peairs

, Choi

, Stewart

, et al. Screening for breast cancer. Semin Oncol, 2017; 44:60.

Sun

, Zhao

, Yang

, et al. Risk factors and preventions of breast cancer. Int J Biol Sci, 2017; 13:1387.

DeSantis

, Fedewa

, Goding Sauer

, et al. Breast cancer statistics, 2015: Convergence of incidence rates between black and white women. CA Cancer J Clin, 2016; 66:31.

Tao

, Shi

, Lu

, et al. Breast Cancer: Epidemiology and etiology. Cell Biochem Biophys, 2015; 72:333.

Eigeliene

, Saarenheimo

, Jekunen

. Potential of liquid biopsies for breast cancer screening, diagnosis, and response to treatment. Oncology, 2019; 96:115.

Saif

, Kasmi

, Allali

, et al. Prediction of DNA methylation in the promoter of gene suppressor tumor. Gene, 2018; 651:166.

Xie

, Kagiampakis

, Pan

, et al. DNA methylation patterns separate senescence from transformation potential and indicate cancer risk. Cancer Cell, 2018; 33:309.

Liu

, Sun

, Qin

, et al. CDH1 promoter methylation correlates with decreased gene expression and poor prognosis in patients with breast cancer. Oncol Lett, 2016; 11:2635.

Noemi

, Idoia

, Natalia

, et al. Differential involvement of RASSF2 hypermethylation in breast cancer subtypes and their prognosis. Oncotarget, 2015; 6:23944.

10.

Martín-Sánchez

, Pernaut-Leza

, Mendaza

, et al. Gene promoter hypermethylation is found in sentinel lymph nodes of breast cancer patients, in samples identified as positive by one-step nucleic acid amplification of cytokeratin 19 mRNA. Virchows Arch, 2016; 469:51.

11.

, Zhanwei

, Lenora Wm

, et al. LINC00472 expression is regulated by promoter methylation and associated with disease-free survival in patients with grade 2 breast cancer. Breast Cancer Res Treat, 2015; 154:473.

12.

Wang

, Guo

, Qian

, et al. miR-506 attenuates methylation of lncRNA MEG3 to inhibit migration and invasion of breast cancer cell lines via targeting SP1 and SP3. Cancer Cell Int, 2018; 18:171.

13.

Shieh

, Hu

, Ma

, et al. Breast cancer risk prediction using a clinical risk model and polygenic risk score. Breast Cancer Res Treat, 2016; 159:513.

14.

, Feng

, Miron

, et al. Breast cancer risk prediction using a polygenic risk score in the familial setting: A prospective study from the Breast Cancer Family Registry and kConFab. Genet Med, 2017; 19:30.

15.

Rudolph

, Song

, Brook

, et al. Joint associations of a polygenic risk score and environmental risk factors for breast cancer in the breast cancer association consortium. Int J Epidemiol, 2018; 47:526.

16.

Ritchie

, Phipson

, Wu

, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res, 2015; 43:e47.

17.

Eberly

. Correlation and simple linear regression. Methods Mol Biol, 2007; 404:143.

18.

Wang

, Wang

, Hang

, et al. A novel gene expression-based prognostic scoring system to predict survival in gastric cancer. Oncotarget, 2016; 7:55343.

19.

Goeman

. L1 penalized estimation in the Cox proportional hazards model. Biom J, 2010; 52:70.

20.

Horne

, Garton

. Likelihood cross-validation versus least squares cross-validation for choosing the smoothing parameter in kernel home-range analysis. J Wildlife Manage, 2011; 70:641.

21.

Eng

, Schiller

, Morrell

. On representing the prognostic value of continuous gene expression biomarkers with the restricted mean survival curve. Oncotarget, 2015; 6:36308.

22.

Robin

, Turck

, Hainard

, et al. pROC: An open-source package for R and S+ to analyze and compare ROC curves. BMC Bioinformatics, 2011; 12:1.

23.

Schröder

, Culhane

, Quackenbush

, et al. survcomp: An R/Bioconductor package for performance assessment and comparison of survival models. Bioinformatics, 2011; 27:3206.

24.

Kohl

, Wiese

, Warscheid

. Cytoscape: Software for visualization and analysis of biological networks. Methods Mol Biol, 2011; 696:291.

25.

Wei

, Sherman

, Lempicki

. Bioinformatics enrichment tools: Paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids Res, 2009; 37:1.

26.

Cimpean

, Tamma

, Ruggieri

, et al. Mast cells in breast cancer angiogenesis. Crit Rev Oncol Hematol, 2017; 115:23.

27.

Janvilisri

, Leelawat

, Roytrakul

, et al. Novel serum biomarkers to differentiate cholangiocarcinoma from benign biliary tract diseases using a proteomic approach. Dis Markers, 2015; 2015:105358.

28.

Hua

, Zhong

, Xia

, et al. RBM24 suppresses cancer progression by upregulating miR-25 to target MALAT1 in nasopharyngeal carcinoma. Cell Death Dis, 2016; 7:e2352.

29.

Chen

, Liu

, Zhao

, et al. Induction, modulation and potential targets of miR-210 in pancreatic cancer cells. Hepatobiliary Pancreat Dis Int, 2012; 11:319.

30.

Zhao

, Guo

, Sun

, et al. Methylcap-seq reveals novel DNA methylation markers for the diagnosis and recurrence prediction of bladder cancer in a Chinese population. PLoS One, 2012; 7:e35175.

31.

Liu

, Han

, Shi

, et al. TUFT1 promotes metastasis and chemoresistance in triple negative breast cancer through the TUFT1/Rab5/Rac1 pathway. Cancer Cell Int, 2019; 19:242.

32.

Liu

, Chen

, Sun

, et al. TUFT1 Promotes Triple Negative Breast Cancer Metastasis, Stemness, and Chemoresistance by Up-Regulating the Rac1/beta-Catenin Pathway. Front Oncol, 2019; 9:617.

33.

Vad-Nielsen

, Jakobsen

, Daugaard

, et al. Regulatory dissection of the CBX5 and hnRNPA1 bi-directional promoter in human breast cancer cells reveals novel transcript variants differentially associated with HP1alpha down-regulation in metastatic cells. BMC Cancer, 2016; 16:32.

34.

Liang

, Lin

, Chen

, et al. Prognostic values of distinct CBX family members in breast cancer. Oncotarget, 2017; 8:92375.

35.

Lindqvist

, Farkas

, Wingren

, et al. DNA methylation pattern of the SLC25A43 gene in breast cancer. Epigenetics, 2012; 7:300.

36.

, Guo

, Wang

, et al. TOX high mobility group box family member 3 rs3803662 and breast cancer risk: A meta-analysis. J Cancer Res Ther, 2018; 14(Suppl):S208.

37.

Ding

, Qi

, Shang

, et al. Docking of CDK1 with antibiotic drugs revealed novel therapeutic value in breast ductal cancer in situ. Oncotarget, 2017; 8:61998.

38.

, Huang

, Ma

, et al. GSTT1 and GSTM1 polymorphisms predict treatment outcome for breast cancer: A systematic review and meta-analysis. Tumour Biol, 2016; 37:151.

39.

Martin-Sanchez

, Mendaza

, Ulazia-Garmendia

, et al. CHL1 hypermethylation as a potential biomarker of poor prognosis in breast cancer. Oncotarget, 2017; 8:15789.