DKKs transcripts with molecular subtypes (PAM50) in patients with breast cancer were examined using the GOBO platform, which was developed by Ringnér et al at Lund University Canceromics Branch (Lund, Sweden), which consists of three main modules, including gene expression analysis, gene co-expression and sample prediction. It is a publicly available platform, including breast cancer subtype-related gene expression data and corresponding clinical data from 1,881 subjects with breast cancer based on Affymetrix microarray analysis (http://co.bmc.lu.se/gobo). The breast tumor set from the 1,881 samples comprises 11 public datasets, including Basal (n=280), HER2⁺ (n=214), Lum A (n=410), Lum B (n=414) and Normal-like (n=255) subtypes according to the breast cancer PAM50 molecular subtype classification (15). The expression analyses of DKKs transcripts in the subtypes and histological grade were performed using one-way analysis of variance (one-way ANOVA). Besides, the DKKs associated with the estrogen receptor (ER)-status were analyzed using the Student's t-test. P-values <0.05 were considered to indicate statistically significant differences.

The transcript levels of DKKs associated with breast cancer outcome, together with overall survival (OS) and relapse-free survival (RFS), were analyzed accurately and in detail using the Kaplan-Meier Plotter online platform, which integrates a number of tumor-related microarrays and detailed clinical prognostic information, as previously described (16). DKKs transcripts with molecular subtypes and clinical outcome from patients with breast cancer were analyzed using the GOBO outcome analysis module. The association of clinical outcome with DKK2 gene expression levels in the subgroups of breast cancer was analyzed through OS and RFS as the endpoint and 10-year censoring. Breast cancer datasets were stratified into 3 quantiles based on DKK2 expression (lower, medium and upper quartile), DKK2_low (log2 expression −5.749 to −0.305), DKK2_ medium (−0.305 to 0.502), and DKK2_high (0.502 to 5.422). Subsequently, the Cox proportional hazards model was applied to estimate the hazard ratio, and the ER status, node status, grade, age and tumor size indicator were used as covariates. Logrank P-values are shown as −log10 (P-value). P-values <0.05 were considered to indicate statistically significant differences.

We also used our data to further confirm the results of bioinformatics analysis on DKKs and the subtypes. Tissue specimens (17 Normal-like subtype and 28 ER-positive subtype) were obtained from 45 patients (no chemotherapy or radiotherapy prior to resection) at Shengjing Hospital of China Medical University (Shenyang, China) from January, 2011 to September, 2013. Informed written consent was obtained from all participants and the study was approved by Ethics Committee of Shengjing Hospital. Total RNA was extracted from these specimens using the RNeasy Mini kit (Qiagen China, Shanghai, China) according to the manufacturer's instructions. RT-qPCR was carried out using the ABI 7500 Real-Time PCR system (Applied Biosystems, Foster City, CA, USA). The RevertAid First Strand cDNA Synthesis kit (Thermo Fisher Scientific, Burlington, ON, Canada) and SYBR^® Premix ExTaq™ II Kit (Takara Biomedical Technology, Dalian, China) were used. Oligonucleotide primers for DKK2 were used as follows: Forward, 5′-CCCCGTTCATTCCTGTTTG-3′ and reverse, 5′-TTCTCCACGGTCCAATCCT-3′. The primers for the internal control, GAPDH, were forward, 5′-CAATGACCCCT TCATTGACC-3′ and reverse, 5′-TGGAAGATGGTGATGGG ATT-3′. The amplification of DKK2 and GAPDH was performed with 1 cycle at 95°C for 10 min and 40 cycles of 95°C for 15 sec and 60°C for 60 sec. Relative mRNA abundance was normalized to the internal standard, GAPDH, using the 2^−∆∆Cq method (17). The samples were divided into 2 groups at the DKK2 median and classified as high or low DKK2.

The TCGA project, which was developed for cancer research, integrates gene expression profiles, methylation, copy number variation and mutation information in 33 cancer types, which helps us to understand multi-dimensional genetic alterations in tumors (18). TCGA breast cancer dataset (19) is available at the TCGA data portal, such as cBioPortal for Cancer Genomics (http://www.cbioportal.org/) (20). According to the PAM50 classification, TCGA breast cancer dataset (19) was divided into 5 subtypes, including the Basal (n=98), HER2⁺ (n=58), Lum A (n=230), Lum B (n=125) and Normal-like (n=8) subtypes. Moreover, METABRIC breast cancer dataset covers several gene expression cohorts, and is widely used for the investigation of breast cancer and subtypes. METABRIC breast cancer dataset (21,22) was grouped into the Basal (n=199), HER2⁺ (n=220), Lum A (n=679), Lum B (n=461) and Normal-like (n=140) subtypes based on the PAM50 classification.

As mentioned above, the GOBO database consists of 11 datasets, which contribute to instability in performance and efficiency of the result. Therefore, we used the publicly available breast cancer dataset with homogeneity and a large sample size to verify the results of GOBO. Subsequently, we undertook breast cancer from TCGA (19) and breast cancer from METABRIC (21,22) to validate expression profile of DKK2 in breast cancer subtypes (22). The expression of DKK2 in the PAM50 subtypes was examined using one-way ANOVA with Dunnett's test with GraphPad prism V7.0 software. P-values <0.05 were considered to indicate statistically significant differences.

We analyzed the survival outcomes of breast cancer and the subtypes according to the differential expression levels of DKK2. Kaplan-Meier survival curves with hazard ratio and logrank P-value were calculated and plotted with the online platform, which is available at Kaplan Meier plotter (http://kmplot.com/analyisis). Kaplan-Meier for breast cancer assessed the relevance of the expression levels of various genes with regards to the clinical outcome (23). Furthermore, we examined METABRIC breast cancer samples (21,22) to validate the clinical outcome of DKK2 in the Normal-like subtype. Data were exported from the cBioPortal platform, and GraphPad prism V7.0 software was then used to perform the survival analysis. The breast cancer sample sets were stratified into 2 quantiles based on DKK2 expression as follows DDK2_medium (medium, 0.424), DKK2_low (−3.04 to 0.424) and DKK2_high (0.424 to 0.534).

To clearly elucidate the underlying mechanisms through which DKK2 is involved in the occurrence and progression of breast cancer at the protein level, we selected breast cancer-related proteins by means of known literature retrieval. Specifically, we conducted the following retrieval expression: 'breast cancer^*' and 'gene^*' or 'protein^*'. Subsequently, we performed PPI network analysis via the Search Tool for the Retrieval of Interacting Genes (STRING v10.5, https://string-db.org/), which provides protein interaction from multiple aspects of evidence, including text-mining, experiments, database, gene co-expression, gene neighborhood, gene-fusion and literature co-occurrence (24). The necessary parameters were set as follows: Meaning of network edges, confidence; minimum required interaction score, medium confidence (0.400). After the PPI network was constructed, Cytoscape software (v3.1.6) was used to visualize the network (25).

Gene ontology (GO) and pathway analysis were performed using the functional enrichment analysis tool (FunRich v3.1.3) software (http://www.funrich.org/), which classifies the biological functions, as well as the signaling pathway of the various proteins involved in the protein networks via hypergeometric distribution algorithm (26). GO analysis consists of molecular function (GO_MF), biological process (GO_BP) and cellular component (GO_CC) (27). Pathway analysis was integrated from the Kyoto Encyclopedia of Genes and Genomes database (KEGG). A threshold value of 0.05 was established for P-values and false discovery rate (FDR) (P<0.05, FDR <0.05).

DKK2 mutation analysis was carried out using the Catalog of Somatic Mutations in Cancer database (COSMIC) (https://cancer.sanger.ac.uk/cosmic). The pie charts present detailed information about distribution of mutations types and substitutions on the coding strand in breast cancer. The alteration frequency of DKK2 mRNA in breast cancer was performed using BioPortal. All operations were performed according to the cBioPortal instructions. The database query was based on mutation and altered expression of the DKK2 in TCGA breast cancer dataset (19) and METABRIC breast cancer dataset (21,22).

The analysis of DKK2 mRNA expression with DNA methylation and copy number was conducted using the cBioPortal database. DKK2 copy number was analyzed by GISTIC, which is the tool to identify gene targeted by somatic copy-number alteration. Moreover, DKK2 methylation was analyzed by Human Methylation 450 (HM450) BeadChip kit. cBioPortal has integrated DNA methylation analysis and copy number analysis in a user-friendly manner.

The statistical significance of expression of DKKs transcripts in the subtypes was examined by one-way ANOVA. In statistics, Dunnett's test is a multiple comparison method developed by Canadian statistician Charles Dunnett to compare each of a number of treatments with a single control. Multiple comparisons to a control are also referred to as many-to-one comparisons. After ANOVA, Dunnett's test, a post hoc test, was used to perform multiple comparison test. The association of the DKKs with ER-status and grade was analyzed by the Student's t-test. Survival curves were plotted using Graphpad Prism software v7.0 (GraphPad Software Inc., San Diego, CA, USA), according to the method of Kaplan and Meier, and the (two-sided) log-rank test was used to assess the statistical significance of differences in OS and RFS among distinct groups of patients. Multivariate analysis of prognostic factors for OS and RFS was performed using the Cox stepwise regression model. P-values <0.05 and FDR <0.05 were considered to indicate statistically significant differences.

Gene expression analysis derived from the GOBO expression module revealed that DKK1 expression was significantly higher in the HER2-enriched, Basal and Normal-like subtype, while it was downregulated in the Lum A and Lum B subtypes (P<0.00001, Fig. 1A). The patients with an ER-negative status exhibited a higher DKK1 expression than those with an ER-positive status (P<0.00001, Fig. 1B). The highest DKK1 transcript expression was observed in Grade 3 tumors, compared with Grade 1 and Grade 2 tumors (P<0.00001, Fig. 1C). The expression of DKK2 and DKK3 was significantly higher in the Lum A and Normal-like subtypes (P<0.00001, Fig. 1A). However, the downregulation of DKK2 and DKK3 was observed in the Basal, HER2⁺ and Lum B subtypes. Unlike DKK1, the patients with an ER-positive status exhibited a higher DKK2 and DKK3 expression than those with an ER-negative status (P<0.00001, Fig. 1B). DKK4 expression was higher in the normal-like subtype and lower in the Lum A subtype (P<0.00001, Fig. 1A). A higher expression level of DKK2, DKK3 and DKK4 was observed in the lower grade and ER-positive tumors (P<0.00001, Fig. 1B). These observations provide further evidence of the DKK1-4 mRNA levels associated with the subtypes, ER status and Grade in breast cancer. Notably, we noted a consistent increase in the expression of DKK family members, with the exception of DKK1, in the normal-like subtypes.

The GOBO outcome module suggested that the DKK2 level was associated with lymph node (LN)-negative (P<0.01) status, ER-positive/LN-negative status (P<0.01), Grade 2 (P<0.01) and Normal-like subtype (P<0.001) (Fig. 2). The mRNA expression levels of DKK1, DKK3 and DKK4 were not associated with the clinical prognosis of the breast cancer subtypes (PAM50 classification) and the results were statistically significant. A high DKK2 expression was associated with a good OS of patients with all subtypes of breast cancer (P=0.0477, Fig. 3A). Moreover, a low DKK2 expression was also associated with a poor RFS of patients with all subtypes of breast cancer (P=0.00493, Fig. 3B). On the contrary, a low DKK2 expression predicted a good OS for patients with the normal-like subtype (P=0.00398, Fig. 3C). In addition, a low DKK2 expression predicted a good RFS for patients with the normal-like subtype (P=0.00677, Fig. 3D). On the whole, we comprehensively analyzed the expression profiles and prognostic value of DKKs in breast cancer subtypes. The main data indicated that a high DKK2 mRNA expression was associated with a better prognosis in breast cancer, which is consistent with the findings of a previous study (28). Nevertheless, we focused our attention on the clinical prognostic value of DKK2 in the normal-like subtype.

To verify the accuracy of the expression profiles of DKK2 in the PAM50-based molecular derived from GOBO in subtypes, the TCGA breast cancer dataset (19) and METABRIC breast cancer dataset (21,22) were used to validate the expression profiles of DKK2 in breast cancer subtypes. The METABRIC breast cancer datasets also indicated that the Lum A, Lum B and Normal-like subtypes the expression of DKK2 exhibited an increasing trend (Fig. 4A). Moreover, DKK2 expression was downregulated in the basal and HER2⁺ subtypes (P<0.0001; Fig. 4A). The TCGA breast cancer dataset (19) also indicated that the transcript level of DKK2 was upregulated in the Lum A, Lum B and normal-like subtypes, and downregulated in the basal and HER2⁺ subtypes (P<0.0001) (Fig. 4B). These findings are mostly consistent with the results of the GOBO dataset. Moreover, Dunnett's multiple test showed that the expression of DKK2 in Normal-like was significantly higher than Basal and HER2⁺ subtype (P<0.0001) (Fig. 4A and B). However, there are a few samples for the normal-like subtype in the TCGA dataset, which may account for the fact that the results shown in Fig. 4B are slightly inconsistent with the results shown in Fig. 1A. However, the prognostic role of DKK2 in the subtypes of breast cancer still requires further investigations with larger sample sets in the future.

The Kaplan-Meier platform revealed that a high DKK2 expression was a good prognostic factor for OS in the all subtypes (Fig. 4C). Moreover, an increased DKK2 expression was also an excellent factor for RFS (Fig. 4D). METABRIC breast cancer (21,22) indicated that a downregulated DKK2 mRNA expression was associated with a better prognosis than the upregulation of its expression in the Normal-like subtype (P<0.05; Fig. 4F). These findings are in line with the results obtained from the GOBO platform.

RT-qPCR analysis of DKK2 mRNA expression was carried out in the Normal-like subtype and ER-positive subtype of breast cancer. The results revealed that the mRNA level of DKK2 was lower in both types than in the normal tissues (P<0.05, Fig. 4E). Kaplan-Meier analysis revealed that an elevated DKK2 mRNA expression was associated with an unfavorable prognosis of patients with the Normal-like subtype of breast cancer (P<0.05, Fig. 4G). Through our own data validation, our results are consistent with the results obtained from TCGA, METABRIC and GOBO.

Based on the information from the literature retrieval, we obtained 39 breast cancer-related proteins, including TP53, ESR1, PGR, ERBB2, CTNNB1, AKT1, APC, ARID1A, ARID1B, ARID2, ASXL1, BAP1, BRCA1, BRCA2, CASP8, CDH1, CDKN1B, CDKN2A, CCND1, MDM2, ZNF217, MYC, FGFR1, ZNF703, GATA3, KRAS, MAP2K4, MAP3K1, MAP3K13, KMT2C, NCOR1, NF1, PIK3CA, PTEN, SETD2, SMAD4, SMARCD1, STK11 and TBX3. Furthermore, the PPI network consisted of 40 nodes and 76 edges (average node degree of 1.15 and average local clustering co-efficient of 0.313) (Fig. 5A). We were aware that DKK2 directly interacted with AXIN1, KRAS, MYC, TP53, CCND1, CDH1, CDKN2A and CTNNB1. The strength of the connection between DKK2 and CTNNB1 was the strongest based on the STRING database. Among these genes, TP53, APC and PTEN are typical tumor suppressor genes (29-31). From the analytical results of PPI network, we also found that DKK2 interacted with SMAD4, MAP3K1 and MAP3K11, which referred to TGF-β/Smad4 and MAPK signaling pathways through CTNNB1. Additionally, high frequency mutations in the tumor suppressor mentioned genes above have been detected in breast cancer (32). These findings suggest that DKK2 may play a role in multiple signaling pathways. The tumor suppressive effect of DKK2 depends on the synergistic effect of the tumor suppressor genes.

To further explore the protein functions derived from the PPI network, FunRich was used to analyze functional and pathway enrichment. The top significant terms of GO enrichment analysis in FunRich were shown in context. Among these biological process, protein metabolism, transcription, apoptosis, regulation of gene expression, cell communication, signal transduction and regulation of nucleobase. Besides, lipid kinase, cell adhesion, kinase binding, and transcript factor activity were the major molecular function. These proteins were related to several cellular components, including protein complex, cytoplasm, cytosol, nucleoplasm and nucleus (Fig. 5B–D). The top 10 significant terms of pathway enrichment analysis in FunRich were mTOR, PDGFR-β, EGF-receptor, Class I PI3K, Arf6 trafficking, LKB1, regulation of CDC42 activity, CDC42, integrin-linked kinase and AP-1 transcription factor network (Fig. 5E). The detailed results are not shown in the text and are available at https://github.com/shaoyoucheng.

The pie chart illustrated the information of mutations of substitution missense, nonsense, synonymous, insertion and deletion, which were generated by COSMIC. The substitution missense rate was 75% and the substitution synonymous rate was 25% of the mutant samples of breast cancer. G>A and C>A mutation occupied 100% of the substitution mutation (Fig. 6A and B). The alteration frequency of the DKK2 mutation in breast cancer was analyzed by BioPortal. Less than 0.1% of the mutations in the patients with breast cancer were screened (Fig. 6C). Moreover, the event that genetic mutation for DKK2 results in nonsense alteration in the 230th amino acid site was detected in breast cancer (Fig. 6D).

The scatter plot displayed that DKK2 methylation was not associated with the DKK2 mRNA level in breast cancer by the cBioPortal database (Fig. 6E and F). In addition, the boxplot provided an overview of the DKK2 mRNA level associated with copy number of deep deletion, shallow deletion, diploid, gain, and amplification, which was generated by the cBioPortal database. Diploid was the cause for the decline in the DKK2 mRNA level in breast cancer (Fig. 6G). On the whole, DKK2 copy number alterations, particularly diploid, account for the elevated DKK2 level in breast cancer.