Development of a circHIPK3-based ceRNA network and identification of mRNA signature in breast cancer patients harboring BRCA mutation

Data collection and pre-processing

The clinical follow-up and expression data of breast cancer patients with BRAC1/2 susceptibility genes BRCA were collected from TCGA database (https://portal.gdc.cancer.gov/), which included RNA-Seq expression profiles and single nucleotide variants (SNVs) processed by mutect2. Samples living longer than 30 days were retained after eliminating those without survival time or follow-up information. We converted ensemble gene IDs to gene symbol IDs, and the gene with multiple gene symbols was expressed as median value.

We downloaded GSE173766 dataset containing circRNA and mRNA sequencing data from GEO (https://www.ncbi.nlm.nih.gov/geo/) database. Moreover, we obtained the gene expression profiles of BRCA mutation samples in GSE103091, GSE16446, GSE20685, GSE20713, GSE42568, and GSE48390 from GEO database. RMA function in “affy” package (Gautier et al., 2004) was used to process each GSE dataset, followed by removal of batch effects using removeBatchEffect function embedded in “limma” package (Ritchie et al., 2015) was utilized to remove. We converted the probe to gene symbol, and removed normal tissue samples. Patients lacking overall survival (OS) and status, follow-up information were excluded, while samples survived longer than 30 days were kept. The batch effect was presented using t-distributed stochastic neighbor embedding (TSNE) (Fig. S1) based on the expression profiles of GSE dataset before and after the removal.

Construction of circHIPK3-based ceRNA network

We obtained hsa_circ_0021603 corresponding gene HIPK3 from the circRNA data in the GSE173766 dataset, and then calculated the correlation of hsa_ circ_ 0021603 with mRNA using Pearson. Positive correlated genes were selected under the threshold of Corr >0.8 and p < 0.05. Subsequently, we predicted hsa_ circ_ 0021603 target miRNAs on the circBank website (http://www.circbank.cn/index.html). Meanwhile, we predicted target genes of miRNAs using miRtarbase, mircode, microT, TargetMiner, TargetScan, miRDB, miRanda, miRmap, starbase, PITA, RNA22, and PicTar. We reserved miRNA-mRNA in at least two databases, and these mRNAs were significantly positively correlated with hsa_circ_0021603. A ceRNA network was constructed using mRNAs, circRNA, and miRNAs.

Functional enrichment analysis

GO analysis on biological process (BP), molecular function (MF), and cellular component (CC) categories as well as KEGG analyses using “WebGestaltR” package (Liao et al., 2019) in R were conducted for functional enrichment analysis.

Developing and validating an mRNA prognostic model

In the ceRNA network, 241 mRNAs were used for additional study. We partitioned the GSE dataset at random into Train and Test datasets with a 1:1 ratio. The Train dataset was used to identify prognostic genes with P 0.05 using univariate Cox regression analysis, followed by LASSO regression analysis to refine the model using “glmnet” package (Hastie, Qian & Tay, 2021). Finally, key prognostic genes were determined by stepwise multivariate regression analysis with stepAIC.

The risk score was calculated in accordance with the formula as follows:

Risk score = ∑βi*ExPi, where i indicated the gene expression levels of prognosis related genes, and β represented Cox regression coefficient. Through surv_cutpoint function in “survminer” package (Kassambara et al., 2017), the optimal cutoff was determined and patients in Train dataset were grouped into low-risk and high-risk groups. Kaplan–Meier curves for each risk group were plotted and the significance of differences was defined by log-rank test. Besides, we conducted receiver operating characteristic (ROC) analysis using “timeROC” package (Blanche, 2015) and evaluated the areas under the ROC curve (AUCs) for 1-year, 3-years, 5-years and 7-years. To validate the robustness of predictive risk model, we calculated the risk score and predicted survival probability in TCGA cohort.

Analysis of clinicopathologic features and SNV between risk groups

We analyzed the differences in human epidermal growth factor receptor-2 (Her2), progesterone receptor (PR), N Stage, estrogen receptor (ER), T Stage between high-risk groups using Fisher’s exact test. Differences in TMB in the TCGA dataset were further investigated. In mutation dataset from TCGA dataset, we screened 8,707 genes with mutation frequency > 3 and screened significant somatic mutation between risk groups using Fisher’s exact test with P < 0.05. Association between risk score and TMB was assessed with Pearson correlation analysis. Additionally, risk score distribution in different clinicopathologic features (age, ER, PR, and Her2, T Stage, N Stage, Stage, M stage). Wilcox tests were used for the analysis of differences between two groups, and the Kruskal-Wallis test assessed group differences.

Gene Set Enrichment Analysis (GSEA) analysis between risk groups

Thereafter, we analyzed biological processes between risk groups through GSEA function in “clusterProfiler” package in TCGA dataset and GSE datasets using HALLMARK pathways in h.all.v7.4.symbols.gmt. Candidates were considered as significant enriched pathways under false discovery rate (FDR) < 0.25 and P < 0.05.

Immune characteristics and pathway characteristics differences between risk groups

The scores of 28 immune cells from previous study (Charoentong et al., 2017) were calculated using single sample GSEA (ssGSEA) method between risk groups in GSE datasets. The immune score was evaluated using ESTIMATE algorithm. Additionally, MCP-counter was employed to assess immune score of 10 immune cells. From previous study (Danilova et al., 2019), we obtained several immune checkpoint genes and compared their distribution in risk groups. In order to study the pathways characteristics of risk groups, we scored the pathways in h.all.v7.4.symbols.gmt and human gene signatures in previous published research (Mariathasan et al., 2018) using ssGSEA and analyzed their correlation with risk score in TCGA dataset.

Response to chemotherapy and immunotherapy of risk groups

The tumor immune dysfunction exclusion (TIDE) algorithm (http://tide.dfci.harvard.edu/) evaluates the M2 subtype of tumor-associated macrophages (exclusion), myeloid-derived suppressor cells, tumor-associated fibroblasts that limit the infiltration of T cells in tumors. Furthermore, two different mechanisms of tumor immune escape, including the dysfunction score of tumor-infiltrating cytotoxic T lymphocytes (dysfunction) and rejection score of CTLs by immunosuppressive factors is also studied with TIDE. A high TIDE score denotes a low rate of ICI therapeutic response. Response prediction to immune checkpoint inhibition (ICI) therapy was performed using TIDE. Additionally, we used the “pRRophetic” package (Geeleher, Cox & Huang, 2014) to calculate the half-maximal inhibitory concentration (IC50) in risk groups both in TCGA and GSE datasets.

Performance of predictive model in immunotherapy datasets

IMvigor210 (Balar et al., 2017), GSE135222 (Kim, Choi & Jung, 2020), and GSE78220 (Hugo et al., 2016) are immunotherapy datasets. Risk score in the above three immunotherapy datasets were calculated and then survival chances were predicted using Kaplan–Meier curves with AUCs. Meanwhile, we used Fisher’s exact test to compare the clinical response distribution amongst risk groups, including partial response (PR), complete response (CR), progressive disease (PD), and stable disease (SD). Risk score distribution in different clinical response groups was compared using Kruskal-Wallis test or wilcox.tests.

Construction of nomogram

In TCGA dataset, a decision tree based on stage, age, T stage, risk type, N stage, M stage was produced to determine risk subtypes and we generated Kaplan–Meier curves of risk subtypes. The significant differences were determined by log-rank test. Univariate and multivariate Cox regression analysis were further used to select predictor of prognosis. Clinical outcomes were evaluated by a nomogram constructed using risk type and clinicopathological features (stage and age, T stage, M stage, N stage). The accuracy and reliability of predictive risk model were analyzed with calibration curves and decision curve analysis (DCA).

Cell culture and transient transfection

MCF-10A, MDA-MB-231, MCF-7 cells were commercially bought from ATCC (Manassas, USA). The cells were cultured in Gibco DEME F-12 medium (Thermo Fisher, Scientific, Waltham, MA, USA). The negative control (NC) and circHIPK3 siRNA (Invitrogen, Carlsbad, CA, USA) were transfected into the cells utilizing Lipofectamine 2000 (Invitrogen, Carlsbad, CA, USA). The target sequences for circHIPK3 siRNAs were ACCATATGTTTATCAAACTCAGT (circHIPK3-si).

Quantitative real-time polymerase chain reaction (qRT-PCR)

TRIzol reagent (Sigma-Aldrich, USA) was applied to extract total RNA from breast cancer cells. The HiScript II Q RT SuperMix (Vazyme, Beijing, China) was used to reverse-transcribe 500 ng of RNA into cDNA. The SYBR Green Master Mix (Thermo Fisher, Scientific, Waltham, MA, USA) was used for qRT-PCR in an ABI 7500 Fast System. For a total of 45 cycles, the PCR amplification settings were at 94 °C for 10 min (min), at 94 °C for 10 s (s), and at 60 °C for 45 s. The internal reference was GAPDH. The primer sequences of GADPH were F: 5′-AATGGGCAGCCGTTAGGAAA-3′and R: 5′-GCCCAATACGACCAAATCAGAG-3′. The primer sequences of circHIPK3 were F: 5′-TTCAACATGTCTACAATCTCGGT-3′and R: 5′-ACCATTCACATAGGTCCGT -3′.

Cell viability

The Cell Counting Kit-8 assay (Beyotime, Jiangsu, China) was performed according to the manufacturer’s protocol for cell viability detection. Different treatment groups of cells were cultured at a density of 1 ×10³ cells per well in 96-well plates. CCK-8 solution was applied at indicated time points. After incubating the cells at 37 °C for 2 h, the OD 450 values of each well were detected using a microplate reader.

Transwell assay

Transwell assays for migration and invasion of MDA-MB-231 and MCF-7 were performed. Briefly, in a nutshell, Matrigel (BD Biosciences, Franklin Lakes, NJ, USA) was coated or left uncoated in chambers before cells (5 104) were added (for migration). In the upper layer, serum-free medium was supplied, and in the lower layer, full DMEM medium. Following a 24-hour incubation period, migrating or invading cells were fixed by 4% paraformaldehyde and stained with 0.1% crystalline violet, counted under a microscope.

Statistical analysis

R software (version 4.0.4) was used in all the statistical analyses, with a two-tailed P < 0.05 being considered as statistical significance.

Construction of a circHIPK3-based ceRNA network

With the threshold of Corr > 0.8 and P < 0.05, we screened 1,358 mRNAs that were positively correlated with hsa_ circ_ 0021603. Subsequently, we predicted 45 hsa_ circ_ 0021603 target miRNAs. Through predicting target genes of 45 miRNAs, we reserved miRNA-mRNA in at least two databases, and these mRNAs were significantly positively correlated with hsa_circ_0021603. Finally, a ceRNA network containing 1 circRNA, 40 miRNAs, and 241 mRNAs was constructed (Fig. 1A). Based on functional enrichment analysis of 241 mRNAs, we found that these mRNAs were closely related with breast cancer, ErbB signaling pathway, Ras signaling pathway, and pathways in cancer (Figs. 1B–1E). The findings suggested that circHIPK3 competing with related mRNA for miRNAs, which may play a regulatory role in the occurrence of cancer.

Development and verification of mRNA prognostic model

According to the ceRNA network, 241 mRNAs that were closely related to breast cancer were used for screening mRNAs for prognostic model construction. Through univariate Cox regression analysis, 40 mRNAs associated with prognosis were screened in the Train dataset (Fig. S2A). LASSO Cox regression analysis was employ to further refine the model. Independent variable coefficients that are progressively heading to zero rose as lambda was steadily increased (Fig. S2B). Figure S2C displays the confidence interval under each lambda from 10-fold cross-validation. 21 mRNAs were found when lambda = 0.0182. With stepwise multivariate regression analysis and stepAIC, a total of 11 mRNAs were identified, including four risk mRNAs and seven protective mRNAs (Fig. S2D).

The risk score was calculated of each patient based on 11 mRNAs (Fig. 2A) as follows: riskScore = 0.652*ANGPTL2−0.289*CD1E−0.282*COLEC12−0.582*FAM117B+0.314* JAG2+0.186*L1CAM−0.405*MSL2+0.718*PSME4−0.42*SH2B3−0.31*TTC39C−0.308* ZBTB20. Survival analysis showed high risk patients exhibited poor prognosis compared with low risk patients (P < 0.0001) with 1-year AUC of 0.67, 3-year AUC of 0.75, 5-year AUC of 0.77 and 7-year AUC of 0.77 in the Train dataset (Fig. 2B). Similar results were also found in the Test dataset and GSE datasets, respectively (Figs. 2C–2D). We further validated the robustness of this model in TCGA dataset. Figure 2E revealed that high risk patients had lower OS (P < 0.0001), progression-free interval (PFI), disease-free interval (DFI) (P = 0.0034), disease specific survival (DSS) (P < 0.0001) (P < 0.0001) than that of low risk patients.

GSEA analysis between risk groups

In GSE and TCGA datasets, EPITHELIAL_MESENCHYMAL_TRANSITION, GLYCOLYSIS, HYPOXIA were significantly enriched in high risk group, while ESTROGEN_RESPONSE, BILE_ACID_METABOLISM, FATTY_ACID_METABOLISM, were significantly enriched in low risk group (Figs. 3A–3B).

Analysis of clinicopathologic features and SNV between risk groups

Thereafter, we assessed the distribution of risk score in N stage, T stage, ER, PR, and Her2, stage. Patients with T3 stage, N2 stage, stage, ER negative, PR negative and Her2 positive status had higher risk score (Fig. 4A). As shown in Fig. 4B, we found that high risk patients exhibited higher proportion of T2 stage, N2 stage, positive ER, negative PR, and positive Her2. We further evaluated the differences in genomic changes between risk groups. High risk patients possessed elevated TMB (P < 0.0001) (Fig. 4C) and had a positive correlation between risk score and TMB (P = 4.18e−09, R = 0.191) (Fig. 4D). In mutation dataset from TCGA dataset, we screened 612 significant mutated genes and displayed the top 20 in Fig. 4E. TP53 (33.73%), PIK3CA (33.3%), CDH1 (14.45%) were the most frequently mutated genes.

Immune characteristics between risk groups

Furthermore, we evaluated the immune landscape in high and low risk groups in GSE datasets. We found that effector memory CD8 T cell, Type 1 T helper cell, activated B cell, effector memory CD4 T cell, activated CD4 T cell, eosinophil, activated CD4 T cell, natural killer cell, CD56 bright natural killer cell were significantly increased in low risk patients (Fig. 5A). Besides, low patients had higher adaptive score and innate score (Fig. 5B), StromalScore, ImmuneScore, and ESTIMATEScore were higher than that of high risk patients (Fig. 5C). Figure 5D deciphered that NK cells, cytotoxic lymphocytes and T cells were distinctly infiltrated in tumors of low patients. Several immune checkpoints such as BTLA, CD200, CD200R1, CD274, CD28, ICOS, IDO2, LAG3 and PDCD1LG2 were overexpressed in low risk patients (Fig. 5E).

Relationship of risk score and pathway characteristics

Meanwhile, we compared the distribution of pathways in TCGA dataset and found that there were 29 significant altered pathways in low and high risk groups (Fig. 6A). Association of risk score with 29 significant altered pathways was further analyzed (Fig. 6B). GLYCOLYSIS, EPITHELIAL_MESENCHYMAL_TRANSITION, and HYPOXIA were dramatically positively correlated with the Risk score. Additionally, risk score was in a positive relationship to DDR, cell cycle, DNA replication, homologous recombination, mismatch Repair (Fig. 6C).

Prediction of response to immunotherapy and chemotherapy of risk groups

Through TIDE algorithm, we observed low risk patients exhibited lower TIDE score, indicating a high response to ICI therapy, while high risk patients had higher exclusion and MDSC both in TCGA and GSE datasets (Figs. 6D–6E). Additionally, low risk patient were sensitive to 14 drugs and high risk patients were more sensitive to 29 anti-tumor drugs in TCGA (Fig. 6F). Low risk patient were sensitive to 47 drugs and high risk patients were sensitive to six anti-tumor drugs in GSE (Fig. 6G).

Performance of predictive model in immunotherapy datasets

The performance of risk model was verified in three immunotherapy datasets. High risk patients had unfavorable prognosis in IMvigor210 (P < 0.0001) and GSE78220 (P = 0.0031) than that of low risk patients with good prediction performance (Figs. 7A–7C). The proportion of PD or PD/SD was significantly increased in high risk group. Besides, patients with PD or PD/SD also exhibited higher risk score (Figs. 7D–7E).

Construction of nomogram

Based on the decision tree, we determined six risk subtypes (Fig. 8A). Figure 8B displayed the different survival probability among six risk subtypes, where patients with subtype C3 were high risk patients while patients with subtype C1 and C2 were low risk patients (Figs. 8C–8D). Figures 8E–8F revealed that risk type was the most important factor related to prognosis. Furthermore, we constructed a nomogram and found that risk score had the greatest impact on survival (Fig. 8G). To evaluate the accuracy of this model, calibration curves for 1, 3, and 5-years were generated (Fig. 8H). The observed curve was highly consistent with predicted curve. Figure 8I displayed nomogram and risk score had the most powerful survival prediction ability.

The cell viability, migration and invasion of breast cancer cells were promoted by circHIPK3

We detected the expression of circHIPK3 in MCF-7 cells by qRT-PCR, and we found that the mRNA expression of circHIPK3 was increased in MCF-7 cells in comparison to normal MCF10A cells (Fig. 9A). We subsequently examined the cell viability of MCF-7 and MDA-MB-231 after circHIPK3 inhibition, and the results showed that after circHIPK3 inhibition, the cell viability of these two kinds of breast cancer cell lines viability was significantly suppressed after circHIPK3 inhibition (Figs. 9B, 9C). We then test the invasion and migration of MCF-7 (Figs. 9D–9F) and MDA-MB-23 (Figs. 9G–9I) cells after circHIPK3 inhibition, and we observed that the invasion and migration ability of the two cell lines was significantly diminished after circHIPK3 inhibition.