Downregulation of the enhancer of zeste homolog 1 transcriptional factor predicts poor prognosis of triple-negative breast cancer patients

Breast cancer (BC) tissue microarrays

A total of three tissue microarrays (TMAs), namely HPreD140Su06 (140 tissue cores), HPreD077Su01 (77 tissue cores), and HPreD050Bc01 (50 tissue cores), were purchased from Shanghai Xinchao Biotechnology Co., Ltd. The quality of the TMA samples was strictly controlled using the following criteria: (I) normal breast tissue or BC with histological classification were pathologically diagnosed using conventional hematoxylin and eosin (HE) staining sections; (II) sufficient clinicopathological information of BC patients was provided, including IHC results of ER, PR, HER2, HER2 fluorescence in situ hybridization, and so on; and (III) the patient had not received any anti-tumor therapy before surgery. TMA samples were excluded: (I) BC samples came from male patients; (II) there were insufficient IHC results for classifying the molecular subtypes of BC; (III) specimens were pathologically diagnosed with a special type of BC, such as medullary carcinoma, mucinous carcinoma, and so on; and (IV) fragmentation was found during IHC staining, thus making it difficult to display the protein expression results of EZH1.

Immunohistochemistry (IHC) staining

A two-step IHC staining method was used in this study. Tissue chips were baked at 65 °C overnight, and were dewaxed for 5 min, then Citrate buffer was used for antigen retrieval under high pressure for 2 min and 30 s. Endogenous peroxidase activity was blocked using 3%H₂O₂. After preliminary experiments, it was determined that the best dilution concentration of the primary rabbit polyclonal antibody anti-EZH1 was 1:100 (https://www.abcam.com/ezh1-antibody-ab64850.html). The Zhongshan Golden Bridge PV-6000 was selected as the secondary antibody. Then, a DAB Kit (20×) (DAB-0031) was used for color rendering. When interpreting the IHC staining results, a semi-quantitative score method was used, which was recorded in Table S1. The raw data of in-house IHC result was recorded in Table S2.

mRNA-expression matrices of BC tissue

mRNA datasets of BC tissue and adjacent normal tissue were accessed from ArrayExpress, cBioPortal, Gene Expression Omnibus (GEO), The Cancer Genome Atlas (TCGA), and the Sequence Read Archive. The keywords searched were (breast cancer OR breast carcinoma) AND “Homo sapiens.” The inclusion standards were as follows: (I) the study type should be an expression profiling by array or by high throughput sequencing, and (II) the clinicopathological parameters of each BC patient should be provided for group division. Duplicate samples and treated patients were excluded. The included mRNA matrices were grouped by GEO platforms and merged into platform datasets. The known batch effect produced in the process of dataset integration was eliminated using the ComBat function with an empirical Bayes framework. The NormalizeBetweenArrays function was used to perform data normalization to certify that the expression values had similar distributions across a platform. The batch-effect-eliminated and normalized platform matrices were matched to the clinic-pathological parameters of each individual. Then the mRNA-expression data of TNBC and non-TNBC tissue were abstracted from batch-effect-eliminated and normalized platform matrices for differential expression and co-expression analysis.

Cancer cell line encyclopedia and dependency map

Cancer Cell Line Encyclopedia (CCLE) sequencing data were downloaded to validate the expression status of EZH1/EZH2 in BC cell lines (https://depmap.org/portal/gene/EZH1?tab=characterization; https://depmap.org/portal/gene/EZH2?tab=characterization). We also explored the clustered regularly interspaced short palindromic repeats (CRISPR) screen data of EZH1/EZH2 in BC cell lines using the Dependency Map (DepMap). Gene effect scores (CERES) were used to evaluate the biological functions of EZH1/EZH2 in BC cell lines. Negative CERES may reflect that the growth of BC cancer cells could be inhibited by the knockdown of EZH1 or EZH2.

Survival data from triple-negative breast cancer subtype (TNBC) cohorts

To evaluate the comprehensive prognostic value of EZH1/EZH2, the overall survival (OS), distal metastasis-free survival (DMFS), and prognosis-free survival (PFS) data of the BC patients were gathered from the Kaplan Meier plotter, which included the prognostic data from 48 GEO gene microarrays and two ArrayExpress datasets (Győrffy, 2021). TNBC patients were subsequently filtered to evaluate the prognostic value of EZH1/EZH2 in TNBC patients. TCGA, METABRIC, and GSE25307 were selected to calculate the pooled hazard ratio (HR) of EZH1, all of which contained no less than 100 TNBC patients.

Prognostic molecular signatures in TNBC patients

Kaplan–Meier survival analysis was performed to screen for potential prognosis indicators in the TCGA TNBC cohort. Univariate Cox regression analysis was used to compute HRs from three large cohorts of TNBC patients independently (i.e., TCGA, METABRIC, and GSE25307). The HRs of each screened prognosis indicator were pooled to determine the prognostic molecular signatures in TNBC patients (risk factor: pooled HR > 1, P < 0.05; protective factor: pooled HR < 1, P < 0.05). Risk score of each patient was computed by using the following equations:

$\mathrm{P}\mathrm{o}\mathrm{o}\mathrm{l}\mathrm{e}\mathrm{d}\mathrm{H}\mathrm{R}=e^{\mathrm{c}\mathrm{o}\mathrm{e}\mathrm{f}\mathrm{f}\mathrm{i}\mathrm{c}\mathrm{i}\mathrm{e}\mathrm{n}\mathrm{t}}$ (e stands for natural logarithm) $$\mathrm{R}\mathrm{i}\mathrm{s}\mathrm{k}\mathrm{s}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}={\sum }_{n=i}^{}\left(coefficien{t}_i\ast expression\mathrm{\_}valu{e}_i\right)$$

Nomogram and time-dependent receiver operating characteristic curves were plotted to predict the risk stratification ability of the prognostic-signatures-based Cox model.

Consensus clustering analysis of TNBC patients

Based on the identified prognostic signatures, TNBC patients in the TCGA cohort were subdivided using a method of consensus clustering. Delta area plot and cumulative distribution function (CDF) plot were utilized to select the optimal k value. The prognostic statuses of the different TNBC subtypes were compared.

Single-sample gene-set enrichment analysis

The immune microenvironments of different TNBC subtypes were portrayed by conducting a single-sample gene-set enrichment analysis (ssGSEA). The inputted immune metagenes were the same as previously demonstrated (Charoentong et al., 2017). Pearson’s correlation coefficients (PCCs) between EZH1 expression and immune infiltrations were calculated.

Somatic-mutation analysis of TNBC patients

The mutation landscapes of EZH1 and EZH2 were explored from TCGA and Metabric using cBioPortal. The somatic mutations of different TNBC subtypes were evaluated by reanalyzing Mutation Annotation Format (MAF) files from TCGA. Somatic interactions, especially mutation co-occurrence or exclusiveness, were studied using a paired Fisher’s exact test.

Differential expression analysis

Differently expressed genes (DEGs) were determined by pooling the standardized mean difference (SMD) value using the mRNA-expression data of TNBC patients, as is aforementioned in the “mRNA-expression matrices of BC tissue” section, compared with non-TNBC patients (Upregulated DEGs: SMD > 0, P < 0.05; Downregulated DEGs: SMD < 0, P < 0.05). All gene symbols containing 61,521 entries were used to filter DEGs in TNBC using a server.

Co-expression analysis

EZH1 or EZH2 co-expressed genes (CEGs) were determined by calculating PCCs using the mRNA-expression data of TNBC and non-TNBC tissue, which have been aforementioned in the “mRNA-expression matrices of BC tissue” section. CEGs were summarized from each platform matrix and were counted, which should be ≥3 (Positive CEGs: PCCs ≥ 0.3, P < 0.05; Negative CEGs: PCCs ≤ $-$0.3, P < 0.05).

Chromatin immunoprecipitation sequencing data analysis

Chromatin immunoprecipitation sequencing (ChIP-seq) data of EZH1 and EZH2 transcriptional target were queried and downloaded from the Cistrome database, one of the largest TF databases. Putative transcriptional targets of EZH1 and EZH2 were filtered by setting scores at ≥1. BigWig files were downloaded for peak identification. Integrative Genomics Viewer v2.8.13 software was utilized to visualize the peaks, where HG38 was selected as the human reference genome.

Prospective signaling pathways in TNBC patients

Gene ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and disease ontology (DO) annotations were conducted after intersecting upregulated DEGs, negative CEGs, and the predicted transcriptional targets. Downregulated DEGs, positive CEGs, and predicted transcriptional targets were also intersected and functionally annotated in the same ways. The mapped GO terms were clustered to identify the most enriched entries. Protein-protein interaction network was constructed to determine the hub gene in the pathway.

Potential TF-miRNA-lncRNA regulatory axis in TNBC patients

BC RNA-seq expression data in the format of fragments per kilobase of transcript per million reads mapped (FPKM) were downloaded from TCGA. FPKM data were transformed into transcripts per million (TPM) data using the following algorithm: $$\mathrm{T}\mathrm{P}\mathrm{M}=\mathrm{e}\mathrm{x}\mathrm{p}\left(\mathrm{l}\mathrm{o}\mathrm{g}\left(\mathrm{F}\mathrm{P}\mathrm{K}\mathrm{M}\right)-\mathrm{l}\mathrm{o}\mathrm{g}\left(\mathrm{s}\mathrm{u}\mathrm{m}\left(\mathrm{F}\mathrm{P}\mathrm{K}\mathrm{M}\right)\right)+\mathrm{l}\mathrm{o}\mathrm{g}\left(1\mathrm{e}6\right)\right)$$

lncRNA TPM expression data were extracted from the BC TCGA matrix. Mature miRNA expression data in the format of reads of exon models per million mapped reads (RPM) were prepared for TNBC and non-TNBC specimens in TCGA. Differentially expressed miRNAs (DEmiRNAs) and differentially expressed lncRNAs (DElncRNAs) were calculated from TNBC tissue in comparison with non-TNBC tissue by using log₂(TPM+1) data (log₂|foldchange| > 1, P < 0.05). Putative binding miRNA targets of EZH1 were predicted by miRTarBase and TarBase V.8, which were intersected by DEmiRNAs. The upstream lncRNAs were forecasted by the encyclopedia of RNA interactomes and LncBase Predicted v.2 (score ≥0.8), which were intersected by DElncRNAs.

Drug repurposing by in-silico analysis

The connectivity map (Cmap) was utilized to identify novel therapeutic small molecules for TNBC patients that targeted the most enriched DEGs. The predicted therapeutic agents for TNBC were strictly limited by using the following criteria: (I) the forecasted drugs should be documented in the DrugBank, or (II) there should be experimental evidence to verify the anti-cancer effect of the forecasted drugs. The predicted small molecules were ranked by a normalized connectivity score, and the first-ranked molecule was used to dock with the transcriptional targets of EZH1. The binding energy threshold was $-$5 kcal/mol, and the lower binding energy, the more stable the molecule.

Statistical analysis

When comparing the expression levels of EZH1 or EZH2 between the two groups, an unpaired two-sample t-test was used. The total number (N), mean (M), and standard deviation (SD) of EZH1 and EZH2 expression values were calculated in TNBC and non-TNBC tissue. SMD was selected to appraise the overall expression levels of EZH1 and EZH2, where SMD > 0 or SMD < 0 implied overexpression and low expression, respectively. If the 95% confidence interval of SMD has no overlap with the zero, the SMD value will be significant. For the calculation of summary SMD and HR, a fixed effect model was selected when I² < 50%. Provided that I² was ≥50%, a randomized effect model should be selected. Begg’s test and Egger’s test were utilized to detect publication bias, where insignificant publication bias reflected the stability of the SMD result (Coelho-Júnior, Trichopoulou & Panza, 2021). We calculated the true positive, false positive, false negative, and true negative rates of EZH1 or EZH2 in TNBC compared with non-TNBC tissue, which was subsequently utilized for drawing summary receiver operating characteristic (SROC) curve. Area under the SROC curve (AUC) was selected to determine the discriminatory ability of EZH1 and EZH2 in TNBC and normal breast tissue, as well as intrinsic subtypes. AUC between 0.70–0.90 indicated a moderate discriminatory ability. A log rank test was performed in Kaplan–Meier (KM) survival analysis. P < 0.05 signified significance.

EZH1 was downregulated in TNBC tissue

Previous studies reported controversial expression levels of EZH1 in BC, so the authors explored the expression pattern of EZH1 in TNBC tissue. Herein, a total of 264 tissue samples containing 177 BC and 87 normal breast tissue samples were eligible for evaluating the protein expression levels of EZH1. All samples had HE stain images to prove the pathological diagnosis (Figs. S1A–S1D). The baseline information of the in-house BC patients and the preliminary statistics of the external mRNA-expression matrices are presented in Tables 1 and 2, respectively. Compared with the non-TNBC tissue, both the EZH1 protein and mRNA-expression levels were significantly decreased in the TNBC tissue (Figs. 1A and 1B). Additionally, with an AUC value between 0.70 and 0.90, the ROC curves implied that there was a moderate discriminatory ability of EZH1 in differentiating between TNBC and non-TNBC tissue, which could be demonstrated by both the quantitative results of in-house IHC (Fig. 1C) and external gene microarrays and RNA-seq (Figs. 1D–1F). The IHC staining results of EZH1 in normal breast, TNBC, and non-TNBC tissue are displayed in Figs. 2–4, respectively, where EZH1 antibody was strongly stained in the nucleus of normal breast epithelial cells. Given the significantly downregulated expression level of EZH1 protein in TNBC, the authors subsequently analyzed the downregulation level of EZH1 in BC subgroups, and it was found that EZH1 protein was lowly expressed in TNBC tissue compared with normal breast, HER2-positive, or luminal BC tissue (Fig. S2).

In summary, EZH1 was significantly downregulated in TNBC tissue, with a SMD of $-$0.59 [ $-$0.80, $-$0.37], and no obvious publication bias or heterogeneity was identified (Begg’s test: continuity corrected P = 0.669; Egger’s test: P = 0.224) (Figs. S3A–S3F). Additionally, EZH1 indicated moderate discriminatory ability between TNBC and non-TNBC tissue, as is shown in SROC (AUC = 0.73 [0.69, 0.76]) and diagnostic likelihood ratios (positive ratio: 1.77 [1.49, 2.11]; negative ratio: 0.44 [0.35, 0.56]). In contrast, EZH2 was significantly upregulated in TNBC tissue (SMD > 0), without obvious publication bias or heterogeneity (Figs. S4A–S4C). Additionally, EZH2 also had a moderate discriminatory ability between TNBC and non-TNBC tissue (area under the curve >0.7) (Fig. S4D).

CRISPR screening analysis of EZH1/EZH2 in BC cell lines

Figures. S5A and S5B showed the expression levels of EZH1/EZH2 in a total of 62 different BC cell lines. According to the CCLE, the expression trends of EZH1/EZH2 in BC cell lines were consistent with that in BC tissue. As is opposite to EZH2, EZH1 was significantly downregulated in BC (Fig. 5A) or TNBC (Fig. 5B) cells in comparison with non-cancerous cells. Unexpectedly, the expression discrepancy of EZH1/EZH2 in TNBC cells was insignificant when compared with non-TNBC cells, which may attribute to the obvious heterogeneity in different TNBC cell lines (Fig. 5C). The CRISPR screening CERES of EZH1/EZH2 were negative in different BC cell lines, including TNBC cells, thus indicating that EZH1/EZH2 may be important in the proliferation and survival of BC cells (Fig. 5D).

EZH1 downregulation presaged poor prognosis in TNBC patients

The prognostic value of EZH1 in BC has not been reported, so the authors explored the potential prognostic implication of EZH1. As is shown in Figs. 6A–6C, a lower EZH1 expression was significantly correlated with worse OS (P = 9.5 $e^{-}$⁵), DMFS (P = 2.3 $e^{-}$⁶), and RFS (P = 2.5 $e^{-}$¹⁶) in BC patients (all with HR <1). Moreover, TNBC patients with lower EZH1 expression were likely to exhibit poorer OS (P = 0.1), DMFS (P = 0.054), and PFS (P = 0.0037) (Figs. 6D–6F). In this setting, EZH1 downregulation may presage poor prognosis of TNBC patients.

In contrast, higher EZH2 expression was significantly correlated with worse OS, DMFS, and RFS in BC patients (Figs. S6A–S6C). However, TNBC patients with lower EZH2 expression were likely to exhibit poorer OS, DMFS, and PFS (Figs. S6D–S6F). Therefore, EZH2 upregulation presaged distinct prognosis in BC and TNBC patients.

A prognostic signature model predicted survival probability of TNBC patients

Because TNBC was a heterogeneous malignancy, the authors tried to establish a model to forecast the survival probability of TNBC patients. Based upon KM survival analysis and HR pooling analysis, a total of six potent prognosticators were identified from TNBC patients (Figs. 7A–7F). Among them, AKAP12, GOLGA1, ITGB3, MATN3, and PCDHB5 were the risk factors for TNBC patients (all with HR >1) while FABP7 was a protective factor for TNBC patients (HR <1) (Figs. S7A–S7F). The risk score equation could be expressed as follows, risk score = 0.3577 * AKAP12 + ( $-$0.1508) * FABP7 + 0.1570 * GOLGA1 + 0.2151 * ITGB3 + 0.2390 * MATN3 + 0.1484 * PCDHB5. Surprisingly, higher risk scores were believed to predict worse OS outcomes of TNBC patients (Fig. 7G). Therefore, a Cox risk model was established using these six prognostic signatures (Fig. 7H), which exhibited an outstanding performance in predicting the prognosis of TNBC patients at 3, 5, and 10 years (all with AUC >0.70) (Fig. 7I). Moreover, the six-signature-based model presented a moderate ability to forecast the prognosis of TNBC patients in the GSE25307 training cohort (AUC at 10 years >0.70) (Fig. 7J).

Prognostic signatures facilitated risk stratification of TNBC patients

The prognostic ability of the identified signatures was presented above; the authors subsequently explored the classification ability of such prognostic signatures. Delta area plot, cumulative distribution function graph, and heatmap of consensus matrix all signified that when k = 3, the clustering result was relatively stable; therefore, a total of three clusters were identified from TNBC patients (cluster 1, cluster 2, and cluster 3) (Figs. 8A–8C).

As is shown in Fig. 8D, the tumor immune microenvironments were compared between different TNBC clusters. Intriguingly, when compared with cluster 1, cluster 2–3 had lower infiltration levels of activated CD4⁺ T cell (P = 0.032). Additionally, there was a positive correlation between EZH1 expression and the infiltration levels of natural killer (P = 0.0448), activated B (P = 0.0317), effector memory CD4⁺ T (P = 0.0051), plasmacytoid dendritic cells (P = 0.0031), and memory B cells (P = 0.0003) (all with correlation coefficients >0) (Fig. 8E). Nonetheless, a negative correlation was detected between EZH1 expression values and the infiltration levels of neutrophils, CD56bright natural killer, and type 17 T helper (all with correlation coefficients <0, and P < 0.05). More importantly, subtype 1 TNBC patients displayed significant survival advantage when compared to subtype 2–3 TNBC patients (P = 0.029) (Fig. 8F). In this context, the identified prognostic signatures may be helpful for identifying a malignant TNBC subtype, which was characterized by low activated CD4⁺ T cell infiltration and poor prognosis.

EZH1 and EZH2 somatic mutations displayed mutual exclusivity in breast cancer tissue

The expression status and prognostic values of EZH1 and EZH2 had been elucidated, and the authors next explored their mutation landscapes in BC and TNBC. As is shown in the Fig. 9A, amplification was the most common mutation type for EZH1 and EZH2 in BC tissue. It was noteworthy that there was a tendency of mutual exclusivity between EZH1 and EZH2 mutations, indicating that EZH1 and EZH2 may play a complementary role for their transcriptional repression in BC. The authors further compared the somatic-mutation types of different TNBC subtypes, and noticed that subtype 2–3 TNBC displayed a higher mutation frequency than subtype 1 TNBC. Additionally, TP53 was the most frequently mutated gene in both subtype 2–3 TNBC and subtype 1 TNBC. For subtype 1 TNBC patients, there was significant co-occurrence between TP53 mutations and tumor suppressor gene PTEN mutations (Fig. 9B). However, TP53 mutations co-occurred with PIK3CA mutations in subtype 2–3 TNBC patients (Fig. 9C). Therefore, the genetic mutations of TNBC were heterogeneous, and a high mutation level of subtype 2–3 TNBC may be responsible for its poor prognosis, compared with subtype 1 TNBC.

EZH1 may mediate competing endogenous RNA mechanisms in TNBC

Subsequently, the authors explored the potential EZH1-miRNA-lncRNA regulatory mechanisms in TNBC. We first identified 61 DEmiRNAs and 156 DElncRNAs from TNBC tissue compared with non-TNBC tissue (Figs. 10A and 10B). After intersecting by the predicted miRNA targets of EZH1, hsa-miR-17-5p, hsa-miR-20a-5p, and hsa-miR-942-5p were determined to be potential miRNA targets of EZH1 (Fig. 10C), all of which were significantly upregulated in TNBC. Then the DElncRNAs were intersected by the binding targets of such three miRNAs, and BLACAT1, LINC00839, and XIST were identified as the putative DEmiRNAs targets. Finally, a potential competing endogenous RNA regulatory network of EZH1 in TNBC was established (Fig. 10D), and the binding sites for EZH1 and hsa-miR-17-5p, hsa-miR-20a-5p, and hsa-miR-942-5p were presented in Fig. 10E.

EZH1 downregulation may promote cell cycle signaling pathway

The molecular mechanisms of EZH1 in TNBC were then investigated. This study identified a total of 8,113 TNBC DEGs, containing 4,215 upregulated DEGs and 3,898 downregulated DEGs in TNBC tissue. Since EZH1 participated in transcriptional repression of its target genes, EZH1 downregulation was anticipated to activate its transcriptional targets. Therefore, we intersected the upregulated TNBC DEGs, EZH1 negative CEGs, and putative EZH1 targets, where a total of 789 targets were identified (Fig. 11A). We clustered the GO terms by similarity and found that the upregulated DEGs and transcriptional target genes were obviously aggregated in the cell-cycle-associated biological process and molecular functions (Fig. 11B). In addition, cell cycle was detected to be the most enriched KEGG pathway (Fig. 11C). More interestingly, such upregulated DEGs and transcriptional targets were mapped for BC disease (Fig. 11D). In all, EZH1 downregulation may activate cell cycle signaling pathway, where CHEK1 was identified as the hub gene (Fig. 11E).

EZH2 upregulation may suppress the response to hormone in TNBC tissue

Since EZH2 was involved in the transcriptional repression of its targets, we intersected the putative EZH2 targets, downregulated TNBC differentially expressed genes, and EZH2 negative co-expressed genes, where a total of 770 targets were identified (Fig. S8A). The downregulated EZH2 transcriptional targets were enriched in response to the hormone GO term and pathways in cancer (Fig. S8B). Furthermore, ESR1 (estrogen receptor 1) was identified as a key gene in response to the hormone (Fig. S8C). The transcriptional factor binding sites for EZH2 and ESR1 was predicted (Fig. S8D).

EZH1 may transcriptionally suppress cell-cycle-associated genes

Genes mapped to the cell-cycle KEGG pathway were selected for further analyses. The authors have elaborated that EZH1 and its transcriptional targets exhibited converse expression trend in TNBC. Additionally, according to the Pearson correlation analysis, a negative correlation was found between the mRNA expression levels of EZH1 and CCNA2 (cyclin A2), CCNB1 (cyclin B1), MAD2L1 (mitotic arrest deficient 2 like 1), and PKMYT1 (protein kinase, membrane associated tyrosine/threonine 1) (Table S3). Based on the public ChIP-seq data, the putative transcriptional factor binding sites for EZH1 and CCNA2, CCNB1, MAD2L1, and PKMYT1 were identified (Fig. 11F). In this setting, it was reasonable to propose that EZH1 might transcriptionally repress such cell-cycle-associated genes. When EZH1 was downregulated in TNBC, cell cycle genes, CCNA2, CCNB1, MAD2L1, and PKMYT1, may be activated. More intriguingly, such potential EZH1 targets were negatively correlated to the immune infiltration levels in TNBC tissue (Fig. S9).

Parthenolide was identified as a potential therapeutic agent for treating TNBC patients

Finally, the authors explored the novel therapeutic small molecules for TNBC patients by targeting cell cycle pathway, and the top four small molecules (i.e., parthenolide, emodic acid, cyclosporin-a, and naproxol) for treating TNBC are predicted in Table 3. Parthenolide molecules were used to dock with the transcriptional targets of EZH1 (i.e., CCNA2, CCNB1, MAD2L1, and PKMYT1) (Fig. 12). The binding energy for parthenolide ligand and CCNB1 protein was the lowest, indicating that CCNB1 seems most effective in conjunction with parthenolide. Moreover, it was observed that parthenolide formed hydrogen bond with the 300^th amino acid of CCNB1 protein, leucine, which also implied the stability of the docking result.