Identification of N6-Methyladenosine-Related LncRNAs for Predicting Overall Survival and Clustering of a Potentially Novel Molecular Subtype of Breast Cancer

We aimed to identify a signature comprising N6-methyladenosine (m6A)-related long non-coding RNAs (lncRNAs) and molecular subtypes associated with breast cancer (BRCA). We obtained data of BRCA samples from The Cancer Genome Atlas database. The m6A-related lncRNA prognostic signature (m6A-LPS) included 10 lncRNAs previously identified as prognostic m6A-related lncRNAs and was constructed using integrated bioinformatics analysis and validated. Accordingly, a risk score based on the m6A-LPS signature was established and shown to confirm differences in survival between high-risk and low-risk groups. Three distinct genotypes were identified, whose characteristics included features of the tumor immune microenvironment in each subtype. Our results indicated that patients in Cluster 2 might have a worse prognostic outcome than those in other clusters. The three genotypes and risk subgroups were enriched in different biological processes and pathways, respectively. We then constructed a competing endogenous RNA network based on the prognostic m6A-related lncRNAs. Finally, we validated the expression levels of target lncRNAs in 72 clinical samples. In summary, the m6A-LPS and the potentially novel genotype may provide a theoretical basis for further study of the molecular mechanism of BRCA and may provide novel insights into precision medicine.


INTRODUCTION
Breast cancer (BRCA) is the most common cancer worldwide and a major contributor to cancerrelated death in women. The increasing incidence and therapeutic challenges due to its heterogeneity have made BRCA a global burden (1,2), especially when triple-negative BRCA (TNBRCA) lacks corresponding molecular targets for therapy (3). Patients with poor prognosis require urgent new therapies; molecular subtype and tumor burden need to be considered in subsequent precision therapy concepts (1).
Among multiple modifications of RNA, N6-methyladenosine (m6A) is the most widespread modification of mRNA and involves all steps of RNA metabolism. It can alter the expression of target genes, making them important regulators of tumor processes (4). Regulators of m6A are divided into three types: writers, readers, and erasers. m6A is a reversible and dynamic RNA modification; crosslinking among the three types of m6A regulators is involved in cancer pathogenesis and progression (5). Recently, many studies have reported advances in understanding the underlying mechanisms of m6A modification in oncogenesis and progression of cancer, including BRCA. For example, ectopic expression of ALKBH5 promotes enrichment of BRCA stem cells under hypoxic conditions (6). By suppressing the upregulation of PD-1, CXCR4, and SOX10 mediated by YTHDF2-degradation, FTO inhibits IFNg-induced killing and the response of anti-PD-1 blockade immunotherapy in melanoma cells (7). METTL3 has oncogenic effects by modulating nonsense-mediated mRNA decay of splicing factors and alternative splicing isoform switches in glioblastoma (8). Bioinformatics research has indicated that aberrations of m6A regulators are associated with poor prognosis of BRCA (9).
Long non-coding RNAs (lncRNAs), defined as transcripts longer than 200 nucleotides, do not encode proteins. Aberration of lncRNA expression appears to play major roles in cancer, such as promoting tumorigenesis and metastasis (10); for example, lncRNA BRCART1 promotes BRCA progression by targeting the microRNA (miR)-1303/PTBP3 axis (11); activation of lncRNA TINCR by H3K27 acetylation promotes trastuzumab resistance and epithelial-mesenchymal transition by targeting miR-125b in BRCA (12). m6A modification in non-coding RNAs has been reported to play a critical role in major normal bioprocesses, and many lncRNAs modified by m6A have been discovered (13). However, the effect of the interaction between m6A modification and noncoding RNAs in tumors remains unclear, and the potential mechanisms of m6A modification of lncRNAs that affect BRCA progression have scarcely been studied. LncRNAs regulated by m6A modification have not yet been identified in BRCA, and they may provide effective therapeutic targets to improve patient prognosis.
In this study, we integrated bioinformatics analyses to identify and validate a robust and stable molecular signature for survival prediction in BRCA. We constructed and verified the prognostic significance of a risk-scoring model based on 10 m6A-related lncRNAs. Based on the m6A-related prognostic lncRNAs, we divided BRCA samples from The Cancer Genome Atlas (TCGA) into different genotypes and explored the differences between the clusters. Furthermore, we constructed a competing endogenous RNA (ceRNA) network to forecast the target miRNAs and mRNAs of these m6A-related prognostic lncRNAs and attempted to predict the potential function of these lncRNAs through functional and pathway enrichment analyses of target mRNAs.

Data Source and Processing
From the TCGA website (https://cancergenome.nih.gov/), RNA expression files (Fragments Per Kilobase of transcript per Million mapped reads [FPKM] normalized), and the corresponding clinical data of BRCA samples and normal samples were obtained. To reduce bias in the statistical analysis in our study, BRCA patients with missing overall survival (OS) values or a follow-up time < 30 days and male BRCA patients were excluded. Finally, we acquired data from 1022 BRCA patients from TCGA. Perl5.30.1 software was used to merge the RNA expression values of each sample into a total matrix file. A list of 23 m6A-related genes was collected based on previous publications, including writers (METTL3, METTL14, METTL16, WTAP, VIRMA, ZC3H13, RBM15, and RBM15B), readers (YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, HNRNPC, FMR1, LRPPRC, HNRNPA2B1, IGFBP1, IGFBP2, IGFBP3, and RBMX), and erasers (FTO and ALKBH5). The expression matrices of 23 m6A-related genes were extracted from the total matrix using R4.0.3 software.

Identification of m6A-Related LncRNAs
From the GENCODE website (https://www.gencodegenes.org/ human/release_29.html), Genome annotation files of Genome Reference Consortium Human Build 38 (GRCh38) were downloaded to distinguish lncRNAs from mRNAs in total matrix. m6A-related lncRNAs were screened using the Pearson's correlation analysis (|Pearson R| > 0.5, p < 0.001). The co-expression network was plotted using the "Cytoscape" software (14). Univariate Cox regression analysis was conducted to identify m6A-related lncRNAs associated with OS, and these lncRNAs were identified as m6A-related prognostic lncRNAs.

m6A-Related Risk Prognostic Signature Model Development
All samples were randomly separated into training and testing sets. To construct a m6A-LPS, we implemented the least absolute shrinkage and selection operator (LASSO) Cox regression analysis via "glmnet" R package in the training set, aiming at dimension reduction to select the most important m6A-related lncRNAs (15,16). The tuning parameter was selected via tenfold cross-validation to avoid overfitting of the predicted signal (17). The risk score for each patient was calculated using the following formula: Coef i * x i where Coef i is the LASSO Cox regression coefficient for the target m6A-related lncRNAs and x i is the expression (FPKM) value of each lncRNA. According to the median risk score in the training set, all BRCA patient samples were then divided into high-risk and low-risk groups.

Validation of the Risk-Scoring Model and Characteristics Identification of Risk Subtype
To verify the predictive function and value of the model, the Kaplan-Meier log-rank test, time-dependent receiver operating characteristic (ROC) curve analysis, univariate Cox regression analysis, and multivariate Cox regression analysis were performed to compare survival between the high-risk and lowrisk groups in the training, testing, and total cohorts. In this part, the "survival", "timeROC", and "survivalROC" packages for R were used for analysis, while the "survminer" and "pheatmap" packages for R were used for drawing plots. Next, we implemented principal component analysis (PCA) using R software for both three datasets and using the "ggplot2" package to plot scatter diagrams.
In each subtype with different clinicopathological characteristics, we applied our model to separate the samples into the high-and low-risk groups to conduct stratification analysis which compared survival between the two risk groups in corresponding clinicopathological subtypes. In the total cohort, the expression of the 10 m6A-related lncRNAs in the model was compared between the high-and low-risk groups. The risk score was compared between subgroups of clinicopathological characteristics, genotypes of m6A-related lncRNAs, tumor mutation burden (TMB) conditions, and TP53 gene mutation status using the Wilcoxon test. The TMB for each sample was calculated using the Perl language. The waterfall grams about the mutation status of top 20 mutational genes and m6A regulators were plotted using the "GenVisR" R p a c k a g e ( 1 8 ) . T h e r e s p o n s e o f c l i n i c a l c o m m o n chemotherapeutics were predicted by the "pRRophetic" R package for risk groups (19).

m6A-Related LncRNAs Genotyping Analysis and Characterization
Based on the expression of m6A-related prognostic lncRNAs, a dimensional reduction of lncRNA genotyping was performed by K-means clustering analysis using the "ConsensusClusterPlus" R package (16,20). We divided n samples into k clusters in which each sample belonged to the cluster with the most similar mean. To determine the construct validity of the genotyping clusters, we employed the Kaplan-Meier log-rank test to compare survival between the clusters. The expression of m6A regulators, tumor microenvironment scores, and TMB were compared using the Kruskal test. The correlations between the different clusters and clinicopathological characteristics were also compared.

Immune Characteristics Analysis of m6A-Related LncRNAs
The expression of human leukocyte antigen (HLA) genes was compared between the different genotypes. The data on the proportions of 22 tumor-infiltrating immune cells in each sample were calculated using xCell, TIMER, quanTiseq, MCPcounter, EPIC, CIBERSORT-ABS, and CIBERSORT for TCGA BRCA samples downloaded from the TIMER 2.0 database website (TIMER (shinyapps.io)). We then analyzed the correlation between tumor-infiltrating immune cells and the expression of the m6A-related lncRNAs and risk score using Spearman's correlation analysis. We normalized the gene expression files of each BRCA sample to predict the proportions of 22 tumor-infiltrating immune cells in each sample via the "CIBERSORT" R package and selected significant results for further analysis. The Kruskal test was used to analyze differences in tumor-infiltrating immune cells between the three genotypes. The difference expression levels of immune check point genes between risk groups were analyzed using the "limma" R package.

Functional and Pathway Enrichment Analysis
Differentially expressed genes (DEGs) between clusters were screened using the "limma" R package with |log2(Fold change)| > 0.5 and p < 0.05 and gene ontology (GO) term functional enrichment analysis were implemented. Gene set enrichment analysis (GSEA) was performed using "org.Hs.eg.db" R package (21) to explore the biological characterization of different genotypes and the low-and highrisk groups were analyzed similarly.
The target miRNAs of the 13 m6A-related prognostic lncRNAs were predicted using the miRcod database (http:// www.mircode.org/) using the Prel programming language. Then, the shared target mRNAs of these miRNAs were screened in miRTarBase, miRD8, and TargetScan databases. DEGs were screened (with |log2(Fold change)| > 0.5 and p < 0.05) using the "limma" R package between BRCA and normal breast tissue (22) to determine the intersection between the DEGs and target mRNAs. We used "Cytoscape" to plot the ceRNA network. To predict the molecular mechanism of m6Arelated lncRNAs, all the mRNAs in the intersection were subjected to functional and pathway enrichment analysis, including GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses using R software, and the loop graph was drawn using the Sangerbox tools (http://www.sangerbox. com/tool).

Quantitative Real-Time Polymerase Chain Reaction for Clinical Samples
In total, 72 BRCA and normal breast samples were collected from patients who underwent surgical treatment at the Department of Breast and Thyroid Surgery, the Third Xiangya Hospital of Central South University from 2019 to 2021. Fortyeight of the 72 samples included 24 luminal BRCAs, 22 Her-2 positive BRCAs, and 2 triple-negative BRCAs (TNBC). The 24 normal breast tissue samples were obtained from BRCA patients who underwent mastectomy. All the fresh samples were stored in liquid nitrogen. This study was approved by the Medical Ethics Committee of the Third Xiangya Hospital of Central South University. Sample acquisition and usage were implemented in accordance with approved guidelines. Each patient signed an informed consent form.
To validate the expression levels of m6A-related lncRNAs in clinical samples, we extracted total RNA from BRCA and normal breast tissue using RNA TRIzol reagent (Shanghai YEASEN Co., Ltd). cDNA of each sample was obtained using a reverse transcription kit (Guangzhou Ribobio Co., Ltd) according to the manufacturer's instructions. The LightCycler 480 Real-Time PCR System was used for quantitative real-time polymerase chain reaction (qRT-PCR) analysis. GAPDH mRNA expression was used as an endogenous control, and the 2 -DDCT formulation was used to calculate related lncRNA expression levels. The primer sequences of the corresponding lncRNA involved in our study are shown in Table 1.

Data Preparation and Identification of m6A-Related LncRNAs
Expression data of 1022 BRCA samples and 112 normal samples were downloaded from TCGA, from which we obtained data for a total of 13,162 lncRNA expression profiles and 19,322 proteincoding gene expression profiles ( Table 2). Next, the expression matrices of 23 m6A regulators were extracted from the total expression matrices. Based on the expression of m6A-related genes, 694 lncRNAs were identified as m6A-related lncRNAs ( Figure 1A). Combined with the prognostic information, univariate Cox regression analysis was performed (p < 0.05) to obtain 13 m6A-related prognostic lncRNAs from the 694 lncRNAs ( Figure 1B). Most lncRNAs were co-expressed with m6A writers, and no lncRNAs were corelated with m6A erasers. Compared with normal breast tissue, besides AL138789.1, BRCA tissue had lower expression levels of these lncRNAs ( Figure 1C). To analyze the effect of the 13 lncRNAs on OS of BRCA, we performed the Kaplan-Meier survival analysis and the results showed that higher expression of LRRC8C-DT, AL51319.1, AL136531.1, and COL4A2-AS1 and lower expression of OTUD6B-AS1 were associated with better OS in the overall cohort of BRCA patients ( Figures 1D-H). Other lncRNAs had no significant effect on OS of BRCA.

Validation of the m6A-LPS Model
BRCA patients were divided into low-and high-risk groups based on the median value of the risk scores in the training set. The BRCA patients in the high-risk group had lower OS rates and shorter OS indicated by Kaplan-Meier survival curves (Figures 2A-C). The risk score and survival status of each BRCA patient were represented by risk curves and scatter plots, respectively, in the training, testing, and total cohorts. Compared with the low-risk group, the high-risk group was associated with higher mortality. The risk heatmap showed that the expression of AL138789.1 and OTUD6B-AS1 increased with increasing risk score, whereas the expression of LRRC8C-DT, AC018926.2, AL51319.1, AL021578.1, ZBTB40-IT1, AC004846.2, AL592301.1, and ZNF197-AS1 decreased with increasing risk score ( Figures 2D-F). The training, testing, and total cohort showed similar results.
The ROC curves demonstrated the promising predictive ability of the m6A-LPS model for BRCA patient survival in the training set (1-year AUC=0.713, 3-year AUC=0.659, 5-year AUC=0.598, and 10-year AUC=0.649), and similar trends were observed in the test set and overall cohort ( Figures 2G-I). These results indicated that the m6A-LPS model had a stable and strong predictive ability for OS in BRCA patients.
A diacritical pattern between the low-and high-risk groups was demonstrated by PCA in the three cohorts. There may be a difference in the m6A modification status of lncRNAs between the risk subgroups ( Figures 3A-C).

Identification and Comparison of Gene Expression of the m6A-Related LncRNAs Genotype
To explore the m6A modification patterns for lncRNAs, Kmeans clustering analysis was performed to obtain appropriate grouping according to the expression of 13 m6A-related prognostic lncRNAs. The total cohort was stratified into three groups: Cluster 1 (n = 675), Cluster 2 (n = 50), and Cluster 3 (n = 297) ( Figures 3D-F). The results of correlation analysis between clinical characteristic and genotypes shown in the heatmap which helped to realize the traits of the populational

The Correlation of m6A-LPS Model and Clinicopathological Features
To verify whether there is a correlation between clinicopathological features and risk scores associated with the m6A-LPS model, we performed the Wilcoxon test and Kruskal-Wallis test, and the results showed that BRCA patients with younger age, higher American Joint Committee on Cancer (AJCC) (23) stage level and TMB, lower immune score and tumor purity, mutation type of TP53, wild-type PIK3CA, and Cluster 2 of m6A-related lncRNA genotypic subgroup had higher risk scores ( Figures 4A-H), whereas the mutation status of BRCA1/2 did not correlate with the risk score ( Figure 4J). Compared with other clinicopathological features, the results of the multi-index ROC curve further supported the strong predictive ability of our model (AUC = 0.717) ( Figure 4I). To further evaluate the prognostic value of the m6A-LPS model, stratification analysis was performed to verify whether the model retained its ability to forecast OS in various clinicopathological subgroups. Next, we confirmed that the m6A-LPS model retained its ability to predict OS for BRCA patients with age older than 35 years, luminal BRCA type, Her-2 positive BRCA type, TNBRCA, higher TMB, TP53 mutant type and PIK3CA wildtype, Cluster 1 of m6A-related lncRNAs genotype, as well as both patients with lower and higher AJCC TMN stage ( Figures 5A-J).
To make the Figure in the correct sequence, we moved this text to the front. These results demonstrated that our model is a stable and potentially predictive tool for patients with BRCA.

Assessment of the m6A-LPS Model as an Independent Prognostic Factor for BRCA Patients
Based on the results of univariate and multivariate Cox regression analyses, we confirmed that the m6A-LPS model could be regarded as an independent prognostic factor for BRCA patients. In contrast, TMB and mutations in TP53, PIK3CA, and BRCA 1 cannot serve as independent prognostic factors for BRCA patients. The higher risk score, younger age, and higher AJCC stage level were significantly associated with worse survival in univariate Cox regression. The statistical significance of risk score and AJCC stage level were further verified in multivariate Cox regression; however, age was excluded ( Figures 5K, L). The m6A-LPS model should serve as a reliable independent prognostic factor for OS in BRCA patients.

Somatic Mutation in Breast Cancer and TMB Value Estimation
First, we extracted the variation in each TCGA BRCA sample to obtain the mutation status of each gene in 872 (88.4%) out of 986 samples. The waterfall diagram revealed that missense mutations, frameshift inserts, and frameshift deletions were common, and we integrated data of m6A-related lncRNA genotypic and risk subgroups to display the somatic mutation status of the top 20 higher mutation rate genes in 789 BRCA samples of the intersection (Figure 6Aa). The somatic mutation status of the m6A genes is shown in the waterfall diagram; similarly, ZC3H13, LRPPRC, YTHDF1, and FMR1 had mutation rates > 10% (Figure 6Ab).

Immune-Related Characteristic Analysis
To investigate the immune characteristics mediated by m6Arelated lncRNAs, we found that m6A-related lncRNAs were closely associated with many immune cells, according to the results of the correlation analysis ( Figure 6B); for example, the  expression levels of LRRC8C−DT were positively correlated with M0 and M2 macrophage abundance but negatively correlated with CD8 + T cell abundance, signifying that the increasing infiltration of CD8 + T cells in the BRCA tissue with a lower expression level of LRRC8C−DT than that in the normal breast tissue. The different expression of HLA genes may reflect coincident immune microenvironment characteristics in m6A-related lncRNA genotypes: Cluster 1 had the highest HLA gene expression, Cluster 2 had the lowest HLA gene expression, and Cluster 3 was in the middle ( Figure 6C). The Kruskal test was performed to compare the tumor-infiltrating immune cells between the clusters; we found a higher abundance of naive B cells and CD8 + T cells in Cluster 1, and a higher abundance of M0 macrophages in Cluster 2 ( Figures 6D-F).
Spearman's correlation analysis was employed to analyze the relevance between the risk score of BRCA patients and the abundance of tumor-infiltrating immune cells, which was calculated using seven software programs, and the statistically significant results (p < 0.05) are displayed in the bubble diagram ( Figure 7A). We found that the abundance of CD4 + and CD8 + T cells was negatively correlated with the risk score for most software. The tumor microenvironment score was calculated based on normalized gene expression files of each sample and was compared using the Kruskal test. Based on the results, we found that Cluster 2 had lower immune and stromal scores but a higher tumor purity; in contrast, Cluster 1 had higher immune and stromal scores but a lower tumor purity (Figures 7B-E). The immune checkpoint genes, including CD200, ICOS, PDCD1LG2, and CD276, had higher expression in the high-  risk group, and the expression of CD44 was higher in the lowrisk group ( Figure 7F). The above results indicate that each group has special immune characteristics and confirm that the modification of m6A methylation lncRNAs plays an indispensable role in the formation of the tumor immune microenvironment in BRCA.

Drug Sensitivity Assessment
To determine the appropriate chemotherapy or other agents for BRCA patients, we performed a drug sensitivity analysis to compare the half maximal inhibitory concentration (IC50) between risk groups. We found that, except for lapatinib, the low-risk group was more sensitive to doxorubicin, docetaxel,  capecitabine, cisplatin, gemcitabine, vinorelbine, and palbociclib ( Figures 8A-H). We suspect that the effect of tyrosine kinase inhibitors (TKIs) is superior in high-risk groups.
GSEA GSEA for the m6A-related lncRNA genotypes was performed to exploit the potential biological processes and pathways leading to molecular heterogeneity between m6A-related lncRNA genotyping. The GO results of comparison in pairs revealed molecular heterogeneity between the 3 m6A-related lncRNA genotyping. We found that many translation-related biological processes (BPs) were enriched in Cluster 1, cell cycle-and mitosis-related BPs were enriched in Cluster 2, and muscle fiber and apoptosis-related BPs were enriched in Cluster 3 (Figures 8Ia-c).
We similarly implemented GSEA for the different risk groups to show that KEGG pathways involving DNA replication, cell cycle, drug metabolism, were significantly enriched in the highrisk group. MiRNAs in cancer and RNA transport were enriched in the low-risk group ( Figure 8J). This may provide some insights into the potential function of m6A-related lncRNAs.

Differential Analysis for m6A-Related LncRNA Genotypes and Functional Enrichment Analysis
To further define the differences in biological characteristics between m6A-related lncRNA genotypes, DEGs in the intersection of three genotypes compared in pairs were excluded, resulting in 717 DEGs ( Figure 9A). We discovered that these DEGs were enriched in GO-BPs, such as signaling receptor activator activity, receptor-ligand activity, and some forms of enzyme activity ( Figure 9B). This provided further insight into the differences between these genotypes.

Construction of the ceRNA Network Based on m6A-Related Prognostic LncRNAs and Functional Enrichment Analysis
To explore the potential function of m6A-related prognostic lncRNAs and how they regulate mRNA expression by sponging miRNAs to participate in BRCA progression, we constructed a ceRNA network based on the 13 m6A-related prognostic lncRNAs that were screened by univariate Cox regression analysis. Three of the 13 lncRNAs from the miRcode database had 48 pairs of interactions with the seven interactive miRNAs. All three lncRNAs were co-expressed with the m6A writer RBM15. In total, 538 target mRNAs of the seven miRNAs were identified from the miRTarBase, miRDB, and TargetScan databases simultaneously. Taking the intersection of all the target mRNAs and the DEGs between BRCA and normal breast tissue together, we obtained 62 target mRNAs to construct the ceRNA network. Finally, our ceRNA network included three lncRNAs, seven miRNAs, and 62 mRNAs ( Figure 10A).
Moreover, functional enrichment analysis, including GO term analysis and KEGG pathway analysis of the 62 target mRNAs, was performed. We found that these genes were enriched in chemorepellent activity, protein phosphatase binding, nuclear receptor activity, and protein tyrosine kinase activity (GO Biological Processes), mitogen-activated protein kinase (MAPK) signaling pathway, cellular senescence, miRNAs in cancer, transcriptional dysregulation in cancer, Ras signaling pathway, etc. (KEGG pathway). All these data could guide us in new directions for further studies to determine the potential functions of these m6A-related lncRNAs in BRCA. In addition, KEGG pathway analysis revealed that these target genes were associated with multiple malignant tumors including prostate cancer, renal cell carcinoma, non-small cell lung cancer, endometrial cancer, acute myeloid leukemia, and glioma ( Figures 10B, C).

Validation of the Expression Level of Four m6A-Related Prognostic LncRNA in BRCA Clinical Samples
According to the difference analysis, the expression levels of m6A-related lncRNAs between BRCA and normal breast tissue  (Figures 10Da-d).

DISCUSSION
Recently, m6A modifications have been reported to not only affect cleavage, transport, stability, and degradation processes of non-coding RNAs, including lncRNAs, but they may also regulate biological cell functions by aberrant expression of lncRNAs (13). These processes may also be involved in diseases such as cancer. In the present study, we screened 13 m6A-related prognostic lncRNAs associated with BRCA including AL136531. for the first time. The lncRNA ZBTB40-IT1 was previously validated to modulate osteoporosis GWAS risk SNPs (rs34920465 and rs6426749) and plays a critical role in bone metabolism that suppresses osteogenesis (24). COL4A2-AS1 can serve as a biomarker in BRCA, and its high expression is related to poor prognosis (25). According to previous reports, OTUD6B-AS1 participates in different mechanisms in multiple tumors, such as inhibiting clear cell renal cell carcinoma proliferation via the Wnt/beta-catenin signaling pathway and targeting corresponding miRNAs to act on thyroid carcinoma, bladder cancer, hepatocellular carcinoma, and colorectal carcinoma (24,(26)(27)(28). Furthermore, we identified a novel molecular signature comprising 10 lncRNAs based on 13 m6A-related prognostic lncRNAs through Cox and LASSO regression analyses and validated them. Kaplan-Meier survival curves indicated a significant divergence in patients who were divided into highand low-risk groups. Time-dependent ROC curves demonstrated that our m6A-LPS was optimal in different cohorts. Moreover, we used a series of analyses to reveal that m6A-LPS can be an independent predictive marker and demonstrated the reproducibility and reliability of m6A-LPS for BRCA prognosis. Unfortunately, we could not find an association between our risk model and the mutation of BRCA1/2. The reasons may be that the mutation rate of BRCA1 and BRCA2 in the general population is only 0.1-0.2% and 0. .4% in all BRCA cases, respectively (29,30). The results of the drug sensitivity analysis showed that the high-risk group was sensitive to lapatinib, which means that the effect of TKI is superior in high-risk groups. In addition, GSEA revealed that the low-risk group was significantly enriched in pathways such as the Janus kinase/signal transducer and activator of transcription signaling pathway, MAPK signaling pathway, Notch signaling pathway, and FC epsilon RI signaling pathway. We believe that m6A-LPS plays a crucial role in the molecular mechanisms of oncogenesis, progression, and prognosis of BRCA. These results confirm that m6A-LPS may provide a reliable prognostic marker and a theoretical basis for the mechanism of BRCA. Despite including training and testing sets to ensure the reliability of results, our study has several limitations. Additional in vitro and in vivo experiments are needed to further confirm the interaction between these lncRNAs and m6A-related genes and how these interactions affect the pathological progress in tumors, especially in BRCA.
In addition, we tried to identify the underlying m6A-related genotype of BRCA and obtained three subtypes (Cluster 1, Cluster 2, and Cluster 3) based on the expression of 13 m6Arelated prognostic lncRNAs. To further explore the different biological characteristics between the genotypes and determine the correlation between these features and m6A-related lncRNAs, we performed a series of analyses and obtained abundance results. We found that most m6A regulators displayed a significantly higher expression status in Cluster 2 than in Clusters 1 and 3. In Cluster 2, the expression of OTUD6B −AS1 and AL138789.1 was higher, while that of AL513190.1 and AC004846.2 was lower. These results suggested that these four lncRNAs are closely correlated with m6A modification. We then assessed tumor-related features in the three clusters. We found that each cluster had unique immune characteristics. Cluster 2 showed a lower immune infiltration level except for M0 macrophages, but a higher tumor purity. The difference analysis of HLA genes showed the same trend: Cluster 2 had the lowest expression level of HLA genes. The test revealed that Cluster 2 exhibited a lower survival risk than other clusters, and the Kaplan-Meier log-rank test showed that Cluster 2 had a worse survival outcome. It is worth noting that there were just 50 patients in Cluster 2. But the correlation analysis between clinical characteristics and genotypes indicated that the BRCA patients with high TMB having a higher proportion in Cluster 2, which would reinforce the association results about the BRCA samples in Cluster 2 with a worse prognosis. Maybe we will amplify the sample size in future analyses to enhance the statistical strength and the reliability of the underlying m6A-related genotype. The role of m6A modification in immunity, especially in infiltrating immunocytes of the tumor microenvironment, has been identified and continues to be researched (31). It has been reported that a lower immune score, which represents a worse immune microenvironment, most likely leads to tumor immune escape with a lower survival rate and a higher recurrence rate (32). This is consistent with the results of our study. This indicates that the m6A-related lncRNA genotype identified in this study may be closely correlated with the immune microenvironment of BRCA. We found that many translationrelated BPs were enriched in Cluster 1, cell cycle-and mitosisrelated BPs were enriched in Cluster 2, and muscle fiber and apoptosis-related BPs were enriched in Cluster 3. DEGs of the three genotypes were enriched in GO-BPs, such as signaling receptor activator activity, receptor-ligand activity, and some kinds of enzyme activity. These results may point to a direction for further research. Finally, we identified target miRNAs and mRNAs for three of the 13 m6A-related prognostic lncRNAs and performed enrichment analysis to explore the potential molecular functions and pathways of m6A-related lncRNAs. The target mRNAs were enriched in pathways such as Ras protein signal transduction and regulation of apoptotic signaling pathway, miRNAs in cancer, further confirming that the three lncRNAs play an important role in the mechanism of BRCA. The mechanism of m6A modification of lncRNAs or the interaction between m6A-related genes and lncRNAs is unclear. It has been reported that m6A modifications may modulate the function of lncRNAs by providing a binding site for the m6A reader proteins or by modulating the structure of the local RNA to induce RNA-binding protein entry, and might also regulate the relationship between lncRNAs and specific DNA sites by affecting the RNA-DNA triple helix structure (13). In the current study, the three m6A-related prognostic lncRNAs involved in the ceRNA network were correlated with RBM15. A recent study showed that XIST could regulate the transcriptional silencing of genes by forming the RNA-binding protein 15 (RBM15)/RBM15B-WTAP-METTL3 complex, which recruits the silencing complex. In addition, knocking down METTL3 or RBM15 reduced the level of m6A modifications on specific transcripts and this resulted in the inactivation of the lncRNA X chromosome (33). Therefore, we believe that the m6A regulator RBM15 plays an important role in lncRNA m6A modification, and the mechanism is worthy of further research.
In conclusion, we identified 13 target lncRNAs associated with BRCA survival; among these, 10 lncRNAs were used to build a prognostic model with m6A modification acting as a novel prognostic clinical trait to predict survival outcome for BRCA patients. Further, we discovered three genotypes relevant to m6A modification lncRNAs, which may provide new insights for precision treatment. Finally, we believe that our results can provide a theoretical basis for further research, and we plan to pursue further studies.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data can be found here: https://cancergenome.nih.gov/.
Hospital of Central South University. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
LQ constructed and supervised this study. BD did the project administration and funding acquisition. XZ performed the data analysis, figures plotted, and writing. XZ and JL did the polymerase chain reaction experiments. JL was responsible for the critical reading of the manuscript. XW performed the figure plotting. XRW and HL were responsible for the data acquisition. All authors read and approved the final manuscript.

FUNDING
This work was supported by the Natural Science Foundation of Hunan Province, China (No. 2020JJ4831).