The Construction of Bone Metastasis-Specific Prognostic Model and Co-expressed Network of Alternative Splicing in Breast Cancer

Background Breast cancer (BRCA) ranks among the top most common female malignancies and was regarded as incurable when combined with bone and distant metastasis. Alternative splicing events (ASEs) together with splicing factors (SFs) were considered responsible for the development and progression of tumors. Methods Datasets including RNA sequencing and ASEs of BRCA samples were achieved from TCGA and TCGASpliceSeq databases. Then, a survival model was built including 15 overall-survival-associated splicing events (OS-SEs) by Cox regression and Lasso regression. The co-expressed SFs of each bone-and-distant-metastasis-related OS-SE were discovered by Pearson correlation analysis. Additionally, Gene Set Variation Analysis (GSVA) was performed to identify the downstream mechanisms of the key OS-SEs. Finally, the results were validated in different online platforms. Results A reliable survival model was established (the area under ROC = 0.856), and CIRBP was found co-expressed with FAM110B (R = 0.320, P < 0.001) associated with the fatty acid metabolism pathway. Conclusion Aberrant SF, CIRBP, regulated a specific ASE, exon skip (ES) of FAM110B, during which the fatty acid metabolism pathway played an essential part in tumorigenesis and prognosis of BRCA.


INTRODUCTION
Breast cancer (BRCA) ranks among the top most common female malignancies in the world (Harbeck et al., 2019). Surgery of the primary tumor remains a cornerstone of curative breast cancer treatment (Mclaughlin, 2013). Systemic therapies for primary breast cancer are quite effective, and adjuvant chemotherapy and adjuvant endocrine therapy are able to decrease the mortality of breast cancer by one third (Early Breast Cancer Trialists' Collaborative Group, Davies et al. 2011. Unfortunately, advanced BRCA is a virtually incurable disease, while metastases are the cause of death in almost all patients, the median overall survival of it is 2-3 years, and the common sites of spread are the bone (most frequent cite), the lungs, and the liver (Cardoso et al., 2018). In order to improve the prognosis of BRCA more efficiently, it is an urgent need to explore the pathogenic mechanism in metastasis and prognosis.
Alternative splicing (AS) took place in human genes quite commonly (Wang et al., 2008) and was critically associated with the carcinogenic process  during which splicing factors (SFs) play key roles. Dysregulating the network built by SFs and alternative splicing events (ASEs), the aberrant AS for some genes and somatic mutations of SFs have been reported that they would cause epithelial-mesenchymal transition and might modulate malignant transformation of cells (Sveen et al., 2016;Kouyama et al., 2019;Wu et al., 2019;Xing et al., 2019). Identification of the aberrant regulation network of SFs and ASEs could not only help predict prognostic and metastatic molecular biomarkers but also find out the potential therapeutic targets (Lee and Abdel-Wahab, 2016;Wang et al., 2019). Most previous studies of BRCA concentrated on alteration at the transcriptome level, but no researchers have studied the analysis of the posttranscriptional process so far.
In this study, based on an all-round analysis of the BRCA dataset, we identified overall-survival-related alternative splicing events (OS-SEs) using Pearson correlation analysis and established a prognostic model by Cox regression accordingly. We then detected the metastasis-related ASEs and their corresponding co-expressed, survival-related SFs and pathways of BRCA. Consequently, the prognostic model we constructed was of great significance for the prediction of BRCA prognosis; we also proposed a potential molecular mechanism and therapeutic target of BRCA.

Data Collection and Preprocessing
This study was approved by the Ethics Committee of The First Affiliated Hospital of Zhengzhou University. The RNAseq data of the data used in this study could be downloaded in the Cancer Genome Atlas (TCGA) database 1 . Demographics, tumor information, and follow-up data of all patients were also retrieved from the database. Furthermore, the seven types of ASE data [alternate acceptor (AA); alternate donor (AD); 1 https://portal.gdc.cancer.gov/ alternate promoter (AP); alternate terminator (AT); exon skip (ES); mutually exclusive exons (ME); retained intron (RI)] along with the Present Spliced In (PSI) values were available in the TCGASpliceSeq database 2 (Ryan et al., 2016). In order to have more reliable results, the data was cleaned under the following rules: (1) Samples containing over 25 of the missing PSI values or without follow-up records were excluded, and the rest of the missing values were supplemented utilizing the K-nearest neighbor algorithm (k = 10).
(2) ASEs whose mean value of PSI less than 0.01 or standard deviations value of PSI less than 0.01 were ruled out.

Independent Prognostic Model Construction
Upset plots were presented to display the relationship between the ASEs and genes directly. For further exploring the contribution of ASE on survival, the univariate Cox regression analysis was adopted using the clinical data to calculate the hazard ratio and its p-value of each filtered ASE. A volcano plot was adopted to show the prognosis-related as well as unrelated ASEs. The top 20 OS-SEs of each splicing pattern were screened according to p-value and were presented in Bubble plots. In addition, to prevent the overfitting of the prognostic model, the selected OS-SEs were used to perform Lasso regression, deleting highly correlated OS-SEs. Accordingly, the final significant OS-SEs were confirmed; n denotes the total number of the selection. The mentioned data of OS-SEs and patents' living status were applied to the multivariate Cox regression model. Thus, the regression coefficient of the jth OS-SE was calculated, denoted as β OS−SEj , j = 1, 2, . . . , n, and the PSI value of the ithsample and the jth OS-SE was represented as PSI OS−SEij , j = 1, 2, . . . , n.
Mathematically, the risk score (RS) of each sample can be described in the following equation: In addition, the samples were then classified into the highrisk group as well as the low-risk group on the basis of the median value of the risk score, and the ROC curve was applied to demonstrate model accuracy. At last, to better compare the difference between the two groups, the survival curve of each group was estimated by Kaplan-Meier survival analysis and the risk curve, scatterplot, and expression heat map were demonstrated as well. Consequently, both univariate and multivariate Cox regression analyses were applied to the processed dataset with age, gender, grade, stage, TNM stage, and risk score. Then, factors whose p-value of both univariate and multivariate Cox regression analyses were less than 0.05 were finally considered as independent prognostic factors.

Discoveries of Relationship Between SF, Metastasis, Stage, and OS-SEs
The dataset with 390 SF factors was achieved from the SpliceAid2 database (Piva et al., 2009). The correlation coefficients and p-values between SFs and OS-SEs were calculated by Pearson correlation analysis to discover their relationship. The SFs and OS-SEs were selected if their absolute value of correlation was more than 0.3 and their p < 0.001. Metastasis information and TNM stage were pivotal paraments of breast cancer so the relationship with OS-SEs may be of great worth. The Kruskal-Wallis test and Mann-Whitney-Wilcoxon test were performed and drawn in Beeswarm plots.
Based on the analyses above, Cytoscape (3.7.1) (Shannon et al., 2003) was used to present the interaction of SF, metastasis, stage, and OS-SEs with high risk as red circles, low risk as purple circles, and the related SF as arrows and the lines indicating their relationships.

Downstream Signaling Pathway Analysis
In order to better compare the results of enrichment analysis and conduct co-expression analysis, quantitative results were quite essential. In this case, Gene Set Variation Analysis (GSVA) (Hanzelmann et al., 2013) was adapted to find out the possible KEGG pathways in BRCA. Prognosis-related ones were selected by univariate Cox analysis. Based on the results above, KEGG pathways co-expressed with metastasis-related OS-SEs were discovered by Pearson analysis. In this way, the potential mechanism of the bone metastasis related to BRCA was identified.

Multidimensional Validation
Genes in the core of the protein-protein interaction networks in ties with the KEGG pathway were detected by Pathway Card 3 . The results of the critical biomarkers were validated to reduce bias conclusions in different dimensions. The significant biomarkers above need external validation to reduce the false-positive rate. Then, a multidimensional validation was performed here. PROGgeneV2 (Goswami and Nakshatri, 2014), Gene Expression Profiling Interactive Analysis (GEPIA) , UALCAN (Chandrashekar et al., 2017), Linkedomics (Vasaikar et al., 2018), cBioportal (Cerami et al., 2012), and Kaplan-Meier plotter (Nagy et al., 2018) showed the relationship between key genes and patients' survival. The Human Protein Atlas (Uhlen et al., 2015) and Genotype-Tissue Expression (GTEx) (Consortium, 2015) demonstrated the expression levels of genes and proteins at the tissue level. Oncomine (Rhodes et al., 2004) presented the results of a multi-study meta-analysis of key genes at the transcriptional level. The Cancer Cell Line Encyclopedia (CCLE) (Ghandi et al., 2019) illustrated the co-expression of the key genes at cellular levels, and String (Snel et al., 2000) displayed the network of SFs and OS-SEs and the key members of the potential pathway. More importantly, they validate in different omics levels. CCLE, cBioportal, and Linkedomics validate at the genomics level; cBioportal, CCLE, UALCAN, Linkedomics, GEPIA, PROGgeneV2, and SurvExpress validate at the clinical level. All datasets except String mentioned in this section perform at the transcriptomic level.

Independent Dataset Validation and Spatial Transcriptome Validation
One independent dataset with a metastasis-free survival (MFS) endpoint was used for independent dataset validation (accession number: GSE11121) 4 (Schmidt et al., 2008).
Additionally, the RNA binding proteins are involved in many biological processes such as RNA maturation, transport, localization, and translation. SFs are a special class of RBPs that can mediate alternative splicing. We have speculated that aberrant SF, CIRBP, regulated a specific ASE, exon skip (ES) of FAM110B, during which the fatty acid metabolism pathway might play an essential part in tumorigenesis, bone metastasis, and prognosis of BRCA. However, the regulation mechanism of SF-mediated AS needs to be demonstrated by molecular biological experiments (e.g., RNA Binding Protein Immunoprecipitation (RIP), Luciferase reporter gene assay). Therefore, spatial transcriptome combined with single-cell RNA sequence (scRNA-seq) data of BRCA was used to identify the cell subtype localization of the key genes in the speculative regulatory mechanism 5 . In terms of quality control (QC), genes with count greater than 1 and expressing in at least three single cells were considered for further analysis. Cells with either fewer than 100,000 transcripts or fewer than 1,500 genes were filtered out.
The Seurat method was applied for integrated data analysis (Butler et al., 2018). After the "sctransform" algorithm for normalization, the "vst" method was used to identify variable genes while the "markvariogram" method was performed to find spatial-specific genes. Variable genes were the input as initial features for principal component analysis (PCA) (Butler et al., 2018). Then, the principal components (PCs) with P < 0.05 were filtered by jackstraw analysis and were incorporated into further t-distributed Stochastic Neighbor Embedding (t-SNE), UMAP (Uniform Manifold Approximation and Projection), and Hematoxylin and Eosin (HE) staining slide to identify cell subclusters (resolution = 0.50) (Chung and Storey, 2015). Only the genes with | log2 fold change (FC)| greater than 0.5 and false discovery rate (FDR) value < 0.05 were selected as differentially expressed genes (DEGs) among cell subclusters. The (spatial) feature plots and violin plots were utilized to illustrate the distribution and expression of DEGs, respectively. Additionally, scMatch (Hou et al., 2019), singleR (Aran et al., 2019), and CellMarker (Zhang X. et al., 2019) databases were used as references for defining each cluster. Furthermore, 50 hallmark gene sets were retrieved from the Molecular Signatures Database (MSigDB) v7.1 6 and Gene Set Variation Analysis (GSVA) was performed to absolutely quantify the signaling pathway activity in each single cell (Hanzelmann et al., 2013;Liberzon et al., 2015).

Statistical Analysis
Statistical analyses were performed on R 3.6.1 software (R Foundation for Statistical Computing, Vienna, Austria) 7 with packages of impute, glmnet, UpSetR, survminer, survivalROC, ggplot2, rms, forestplot, beeswarm, and preprocessCore. P < 0.05 were considered statistically significant except for the selection of SFs and OS-SEs where the p < 0.001.

The Overview of the Analysis and Dataset
The schematic view of this integrated analysis is depicted in Figure 1. Supplementary Table S1 summarizes the baseline characteristics of all patients. Patients were not one-to-one congruent with samples. The total amount of primary BRCA patients was 1,097, and the average survival time of the patients was 1,199 days with an average age of 58.5 years. In addition, there were only 1,080 primary samples; 32 of them underwent bone metastasis in the experimental group, and 1,048 samples did not experience bone metastasis in the control group. 7 www.r-project.org There were 5,433 ASEs found in 2,933 parent genes in the samples with 2,728 ESs including 1,801 genes; 1,304 APs including 518 genes; 444 AAs including 396 genes; 414 ADs including 374 genes; 254 RIs including 229 genes; 251 ATs including 117 genes; and 38 MEs including 37 genes. One ASE could take place in different kinds of genes, and one gene could locate in different kinds of splicing patterns (Figure 2A and Supplementary Table S2). 294 of OS-ASEs were associated with survival status. Ranked by the magnitude of the association strength with survival, ES came first, followed by AP, AA, RI, and AD.
Based on the top 20 OS-SEs, we built a prognostic model for BRCA. For avoiding over-fitting, Lasso plot ( Figure 4A) and the Lambda plot ( Figure 4B) were performed. Finally, 15 OS-SEs  The Upset plots of seven types of prognosis-related alternative splicing events and their parent genes. The lower part of each plot describes the permutations of the alternative splicing events; for each circumstance, the alternative splicing event is included if its corresponding place is filled with a red dot. The upper part of the plot represents the number of genes for each circumstance. AA, alternate acceptor; AD, alternate donor; AP, alternate promoter; AT, alternate terminator; ES, exon skip; ME, mutually exclusive exons; RI, retained intron.
were selected to apply to the multivariate Cox regression analysis. The result of the regression was quite reliable with the ROC area under the curve (AUC) of 0.856 ( Figure 4C). The median value of the risk score was 0.97, thus dividing the whole samples into two groups. The survival analysis results ( Figure 4D) showed the survival probability of each group, indicating a vast difference between them, and the Kaplan-Meier curve demonstrates the reliability of the model by a p < 0.001. The risk curve and scatterplot indicated that samples with higher risk scores had a higher risk of mortality (Figures 4E,F)  In order to verify the independence of the risk score, both univariate and multivariate Cox regression analyses were applied to age, gender, stage, TNM stage, and risk score. With both p-values of the risk score in two analyses less than 0.001 and hazard ratios which were 1.053 [95% confidence interval (CI): 1.029-1.077] and 1.042 (95% CI: 1.023-1.062), the risk score proved to be a well-predicting model (Figures 4H,I).
Consequently, SFs can serve as a predictor regarding survival.

The Detection of the Potential Splicing Regulatory Network
There were 185 pathways detected by GSVA, which had potential prognostic predictive abilities. Then, just 26 pathways were selected as the prognosis-related pathways by univariate Cox analysis. Among them, ubiquitin-mediated proteolysis had the highest score of 3.36 (P = 0.027).
The interacting and correlating relationship between SFs and OS-SEs was illustrated in the network (Figure 5A). During the detection of cancer status-related OS-SEs, 27, 23, and 43 OS-SEs were found related to distant metastasis, bone metastasis, and co-expression, respectively, and only 2 of them were collapsed ( Figure 5B). The two OS-SEs in relation with both metastasis and stage were UNKL-33078-AP and FAM110B-83922-ES. The different expressions of these two OS-SEs between specific metastasis and primary tumor were depicted in boxplots (Figures 5C-F). The correlations between KEGG pathways and the two OS-SEs were showcased in the heat map (Figure 6) according to the result of Pearson analysis. Pyrimidine metabolism was the pathway most connected to UNKL (R = −0.385, P < 0.001), and fatty acid metabolism was the pathway most related to FAM110B (R = 0.223, P < 0.001).

External Validation
The Human Protein Atlas (Supplementary Figure S1A) showed that FAM110B had a higher protein level in cells of BRCA than normal ones while UNKL was not detected in either normal or breast cancer. Thus, the study and conclusions were mainly focused on FAM110B and the related SF, CIRBP (R = 0.320, P < 0.001). The fatty acid metabolism, the most related KEGG pathway of FAM110B, revealed that ACAT2, ACAT1, ALOX15B, DHCR7, and ACAA1 were the five genes among the most centered genes (Supplementary Figure S11A). The Human Protein Atlas (Supplementary Figures S1B-F) demonstrated that all of these five genes have a relatively high level of protein expression.

(Supplementary
In the end, the interaction PPI network of FAM110B, CIRBP, ACAT2, ACAT1, ALOX15B, DHCR7, and ACAA1 from String (Supplementary Figure S11B) was shown. Thus, we hypothesized that aberrant CIRBP regulated FAM110B-83922-ES in ties with the tumorigenesis, metastasis, and prognosis of BRCA via fatty acid metabolism, and the schematic diagram of this scientific hypothesis is depicted in Figure 7.
genes were incorporated into the multivariate Cox model and the risk score for each BRCA patient was calculated by the formula of the multivariate model. The risk line and scatterplot illustrated the distribution of the risk score among all patients (Figures 8B,C). The Kaplan-Meier survival curve showed the significantly prognostic value of the risk score ( Figure 8D, P < 0.001). Furthermore, the accuracy and goodness of fit (GOF) of the multivariate Cox regression model was illustrated by the ROC curve (AUC = 0.710) and the residual plot (Figures 8E,F).
Spatial transcriptome combined with scRNA-seq data of BRCA was used to identify the cell subtype localization of the key genes in the speculative regulatory mechanism. Eleven clusters and nine clusters were identified by t-SNE and UMAP, respectively ( Figure 9A). Invasive ductal carcinoma and intraductal carcinoma in situ as well as normal tissue were clearly demonstrated in HE staining slides ( Figure 9A). The feature plots and spatial feature plots were utilized to illustrate the distribution and expression of CIRBP, FAM110B, ACAA1, ACAT1, and DHCR7, showing that these key genes were highly expressed in the invasive ductal carcinoma tissue (cluster 2, 5, 8) (Figures 9B,C). Additionally, the cell-cycle analysis suggested that most cells highly expressing these key genes were in the phase of G2M and S (Figures 9D,E), and those cell-division-related signaling pathways such as HALLMARK_G2M_CHECKPOINT and HALLMARK_E2F_TARGETS were activated while some tumor-suppressor signaling pathways such as HALLMARK_ P53_PATHWAY and HALLMARK_TNFA_SIGNALING_VIA_ NFKB were downregulated in the invasive ductal carcinoma tissue ( Figure 9F).

DISCUSSION
BRCA is women's most common neoplasm, whose 5-year survival rate is less than 40% in some areas of the world (Akram et al., 2017). The incidence and mortality rates of BRCA are still expected to increase in the next 5-10 years (Anastasiadi et al., 2017). BRCA could be spread to sites like bone, lungs, the liver, and axillary lymph nodes (Harbeck et al., 2019); death is more likely to occur due to metastasis (Akram et al., 2017). Therefore, it is quite essential for us to delve into BRCA wishing for an improvement of the diagnosis and treatment. ASEs frequently happened in eukaryote lineages (Bush et al., 2017), and alternatively spliced isoforms are generated in more than ninety-five percent of multi-exon genes in human beings (Pan et al., 2008;Wang et al., 2008). In addition, ASEs along with SFs are considered to have essential impacts on the development and progression of tumors (Srebrow and Kornblihtt, 2006), and AS changes might become independent oncogenic processes (Climente-Gonzalez et al., 2017). Most previous studies of BRCA are focused on alterations at the transcriptome level, but analysis of the posttranscriptional process is largely ignored. One study found that some ASEs in gene NF-YA are associated with BRCA (Dolfini et al., 2019). Another study also focused on some OS-SEs in BRCA and analyzed some genes associated with OS-SEs (Zhang D. et al., 2019). Despite the significance, metastasis-related ASEs and potential therapeutic targets were ignored, and the relationship between SF and ASEs was also not analyzed.
In this study, univariate Cox regression and Lasso regression selected 15 OS-SEs and a prognostic model was established showing that SFs could well predict survival. Two OS-SEs were identified to be associated with tumorigenesis and metastasis, and their related SFs were selected along with corresponding prognostic KEGG pathways. After multidimensional validation, we finally assumed that CIRBP regulated FAM110B-83922-ES of BRCA via fatty acid metabolism, and five proteins involved in fatty acid metabolism were demonstrated to be highly related to survival as well. The prognostic model could be an effectual tool for doctors in the future, and the assumption might provide an idea for finding a particular target to increase the survival rate of BRCA.
CIRBP, cold-inducible RNA-binding protein, belonging to the cold-shock protein family, is expressed in different cell types and was involved in various cancer pathological processes, including cell survival, cell proliferation, and cancer (Zhu et al., 2016;Su et al., 2020). It was reported that CIRBP had different effects in various cancers according to the cell context such as acting as a tumor suppressor in rectal carcinoma, endometrial carcinoma, and ovarian tumor and having protumorigenic influences on melanoma, colorectal cancer, prostate cancer, central-nervous-system-related tumor, liver-pancreas carcinomas, skin squamous cell carcinoma, bladder cancer, and pituitary corticotrope adenoma (Zhu et al., 2016;Liao et al., 2017;Zhong and Huang, 2017;Chen et al., 2018;Colasanti et al., 2018;Lujan et al., 2018;Lin et al., 2019). Some studies found down-expressed CIRBP to be always related with poor-survival cell-cycle process and disease stage in BRCA patients (Joe and Nam, 2016;Rodrigues-Peres et al., 2019). Additionally, a protein microarray experimental study sheds light on CIRBP as an autoantibody target potentially used for early BRCA FIGURE 8 | Independent dataset validation by GSE11121. In GSE11121, univariate Cox regression analysis suggested that CIRBP (P < 0.001), FAM110B (P = 0.013), ACAA1 (P = 0.045), ACAA2 (P = 0.043), ACAT2 (P = 0.037), and ALOX15 (P < 0.001) had a significant association with metastasis-free survival (MFS) of BRCA patients (A). Then, these key genes were incorporated into the multivariate Cox model, and the risk score for each BRCA patient was calculated by the formula of the multivariate model. The risk line and scatterplot illustrated the distribution of risk score among all patients (B,C). The Kaplan-Meier survival curve showed the significantly prognostic value of the risk score (D, P < 0.001). Furthermore, the accuracy and goodness of fit (GOF) of the multivariate Cox regression model were illustrated by the ROC curve (AUC = 0.710) and the residual plot (E,F).
detection (Lacombe et al., 2014). Studies also identified CIRBP as a potential prognosis-related gene which could be used in BRCA prognostic prediction (Shimizu and Nakayama, 2019). Acting as a cold-inducible RNA-binding protein, CIRBP also had a positive influence on HuR protein levels concerned with BRCA formation (Kotta-Loizou et al., 2016). FIGURE 9 | Spatial transcriptome validation. Spatial transcriptome combined with scRNA-seq data of BRCA was used to identify the cell subtype localization of the key genes in the speculative regulatory mechanism. Eleven clusters and nine clusters were identified by t-SNE and UMAP (A). Invasive ductal carcinoma and intraductal carcinoma in situ as well as normal tissue were clearly demonstrated in HE staining slides (A). The feature plots and spatial feature plots were utilized to illustrate the distribution and expression of CIRBP, FAM110B, ACAA1, ACAT1, and DHCR7, showing that these key genes were highly expressed in the invasive ductal carcinoma tissue (clusters 2, 5, 8) (B,C). Additionally, the cell-cycle analysis suggested that most cells highly expressing these key genes were in the phase of G2M and S (D,E), and those cell-division-related signaling pathways such as HALLMARK_G2M_CHECKPOINT and HALLMARK_E2F_TARGETS were activated while some tumor-suppressor signaling pathways such as HALLMARK_P53_PATHWAY and HALLMARK_TNFA_SIGNALING_VIA_NFKB were downregulated in the invasive ductal carcinoma tissue (F).
A study also found the co-regulated expression of SFs in BRCA, indicating that CIRBP expressed much lower in tumor tissue than in metastatic tissue (Koedoot et al., 2019).
FAM110B, family with sequence similarity 110 member B, localized in centrosomes, and the progress in the G1 phase of the cell cycle is affected by the ectopic expression of its protein (Hauge et al., 2007). Much of the researches focused on the relationship between FAM110B and BRCA. FAM110B played a negative role in the translational regulation of BRCA1, which is a tumor suppressor of BRCA (Dacheux et al., 2013). Previous studies have also demonstrated that the signature of FAM110B along with 40 other genes can serve as a survival predictor for BRCA. The signature was associated not only with age and ER status but also with overall survival and distant metastasis-free survival (Yin et al., 2014). What is more, FAM110B expressed significantly lower in young women compared to older women with BRCA, so that it may play quite different roles in the development of BRCA in young women than in older women (Colak et al., 2013).
Until now, there is no evidence about the direct regulation between CIRBP and FAM110B. Furthermore, they are both involved in cell-cycle progression based on the up-to-date reports (Hauge et al., 2007;Zhu et al., 2016;Su et al., 2020), and we found out that they are co-expressed in the BRCA database. To go through the mechanism of CIRBP regulating FAM110B-83922-ES, fatty acid metabolism was detected as the possible pathway connecting them.
Through divergent mechanisms, endogenous fatty acid status could have a great impact on health and disease (Brenna et al., 2018). Consistent with our study, many studies illustrated the significance of fatty acid metabolism in BRCA. Characteristics of different types of BRCA vary in fatty acid metabolism (Yamashita et al., 2017). In addition, fatty acid metabolism was engaged in BRCA progression, and some proteins related to fatty acid transport revealed to enhance the possibility of migration and invasion of BRCA (Yen et al., 2018). Furthermore, evidence demonstrated that fatty acid metabolism was considered a potential target for BRCA therapy (Kinlaw et al., 2016).
In sum, although the results of our study were quite reliable, our study still had some limitations. The result was only based on one single chart of sequencing data and computer analyses where the sample size of the dataset and the accuracy of molecular mechanisms were limited. However, the most important shortcoming is that systematic errors could not be avoided and wet experiments and clinical trials are required to perform to verify the results. The inherent vice of in silico analysis was bias based on different platforms, and our scientific hypothesis only based on bioinformation instead of mechanism exploring. In addition, a small sample size for bone metastasis samples may result in an unbalance of the datasets. In order to better diminish the bias, external multidimensional validation was followed to minimize the negative impact of this limitation, showing the reliability of our findings.

CONCLUSION
In this comparative bioinformatics analysis, we constructed a prognostic model of BRCA. In addition, our further findings suggest that the aberrant splicing factor, CIRBP, regulated an alternative splicing event, the exon skip of FAM110B, during which the fatty acid metabolism pathway might play an essential part in tumorigenesis and prognosis of BRCA. This scientific proposition might provide direct instruction for the following biological experiments.

ETHICS STATEMENT
This study was approved by the Ethics Committee of the First Affiliated Hospital of Zhengzhou University. The patients/participants provided their written informed consent to participate in this study.

ACKNOWLEDGMENTS
We thank the Cancer Genome Atlas (TCGA) team for use of their data.