Research Paper Volume 12, Issue 16 pp 16238—16254

Prognostic values of GPNMB identified by mining TCGA database and STAD microenvironment

Kunhou Yao1,2, , Lunshou Wei1,4, , Junjie Zhang1,2, , Chenyu Wang1,2, , Chaoyang Wang1,2, , Changjiang Qin1,2, , Song Li1,3, ,

  • 1 Henan University, Kaifeng 475000, Henan, P.R. China
  • 2 Department of General Surgery, Huaihe Hospital of Henan University, Kaifeng 475000, Henan, P.R. China
  • 3 Department of Urology, Huaihe Hospital of Henan University, Kaifeng 475000, Henan, P.R. China
  • 4 Department of Gastroenterology, Huaihe Hospital of Henan University, Kaifeng 475000, Henan, P.R. China

Received: February 28, 2020       Accepted: June 19, 2020       Published: August 21, 2020      

https://doi.org/10.18632/aging.103646
How to Cite

Copyright © 2020 Yao et al. This is an open-access article distributed under the terms of the Creative Commons Attribution (CC BY) 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

The survival rate of stomach adenocarcinoma patients with immune and stromal scores and different clinicopathological features obtained from the TCGA datasets was systematically compared. A list of genes that are correlated with stomach adenocarcinoma microenvironment were extracted using the TCGA database to predict the prognosis and survival. In addition, the differentially expressed genes were extracted by comparing the immune and stromal scores of the groups. The protein-protein interaction network, and functional and pathway enrichment analyses of differentially expressed genes were performed. A total of 8 hub genes were selected from the differentially expressed genes to predict the overall survival and disease-free survival rates. GPNMB was selected from the hub genes based on the survival and prognosis analyses. A nomogram was built by including the potential risk factors based on multivariate Cox analysis. Cell function experiments and xenograft tumors were conducted in vivo to further verify the role of GPNMB in tumor progression. The predicted microRNA, miR-30b-3p, might act as upstream negative regulator and binding to 3’ UTR of GPNMB, confirming by fluorescent enzyme reporter gene experiment. In summary, immune-related scores are crucial factors in the malignant progression of stomach adenocarcinoma and GPNMB acts as a potentially useful prognostic factor for stratification and in developing the treatment strategy

Introduction

Stomach adenocarcinoma (STAD) is one of the most common malignant tumors of the digestive tract. In 2018, there were approximately 1,033,701 new cases and approximately 782,685 deaths associated with STAD worldwide, and most of these were in the locally advanced stage at the time of diagnosis [1, 2]. The mechanism of occurrence and development of STAD is still unclear. The current treatments of STAD include surgery, radiation therapy, chemotherapy, etc. However, the incidence of local recurrence or distant metastasis of gastric cancer is 40% to 70% after surgery, and the side effects associated after undergoing radiotherapy and chemotherapy are obvious [3]. Therefore, the prevention and control of STAD has become an urgent public health issue, and it is necessary to explore the underlying mechanism of STAD progression in order to find new therapeutic and diagnostic targets that can improve the survival rate of patients.

With the concept of "survival with tumor", the tumor immune microenvironment has attracted much attention and taken immunosuppression as its core feature, as it plays a decisive role in the production and proliferation of tumor cells [4]. Tumor microenvironment is a unique environment that involves interactions of tumor cells, including tumor stroma, adjacent cells, various immune cells around blood vessels and other immune-related mediators, with the host. Changes in various components of tumor microenvironment demonstrated important effects on tumor growth, invasion, metastasis and tumor immune tolerance [5]. The occurrence and development of STAD is based on the coordinated evolution of cancer cells and tumor immune microenvironment. Surgery, radiotherapy, and chemotherapy can aggravate immunosuppression, promoting the proliferation of residual cancer cells or the recurrence of tumor cells [6]. Therefore, to improve immune surveillance, restrain immunosuppression and either activate or inhibit immune escape, it is imperative to improve the immune microenvironment of gastric cancer, affecting the growth, progression and outcome of gastric cancer. This becomes a new target for the prevention and treatment of gastric cancer [7]. Immunization and stromal cells are the two major non-tumor components in the tumor microenvironment and has great value in the diagnosis and prognosis assessment of tumors [8]. Yoshihara et al. have designed an algorithm called ESTIMATE, which is used for calculating immune and stromal scores for predicting the infiltration of non-tumor cells using gene expression data from the The Cancer Genome Atlas (TCGA) database [8].

In this study, the immune and stromal scores of patients with STAD were calculated, and the scores in STAD with different clinicopathological features from the TCGA datasets were systematically compared. Next, a list of genes that are correlated with microenvironment was extracted, and predicted the prognosis and survival in STAD patients via TCGA database of STAD cohorts and ESTIMATE algorithm-derived immune scores. Glycoprotein non-metastatic melanoma protein B (GPNMB) was selected from the hub genes based on survival and prognosis analyses. Multivariate cox regression analyses were performed and a nomogram was built with potential risk factors based on a multivariate Cox analysis to predict the survival probability. Furthermore, cell function experiments and xenograft tumors in vivo were performed to verify the role of it in tumor progression.

Results

Association of immune scores and stromal scores with STAD

To investigate the prognostic value of immune-related scores (stromal scores and immune scores) in STAD, all the 317 STAD cases were divided into high and low groups based on their immune and stromal scores to find the association between the survival rate and the scores. The results revealed that patients with low stromal scores demonstrated slightly increased overall and disease-free survival rates when compared to high stromal scores (Figure 1, P < 0.05), while the overall survival rate of low immune scores was much higher than that of high immune scores (Figure 1, P < 0.05). Next, the association between the scores and clinical information were analyzed based on the data obtained from the TCGA database. Among them, the stromal scores and immune scores showed positive correlation with stages T and M as well as grade, indicating that the immune-related scores might contribute to STAD progression, while only stromal scores were positively correlated with stage N (Figure 2, P < 0.005). In addition, the immune-related scores were analyzed according to CD274 expression and TP53 mutation states in STAD. As shown in Figure 2, TP53 mutant cases had lower scores when compared to those with wild-type TP53 and CD274 expressed in high score group are higher than that of low group.

Immune scores and stromal scores are associated with Kaplan-Meier survival in STAD. The OS and RFS curves between high and low groups based on immune scores (A, B). The OS and RFS curves between high and low groups based on immune scores (C, D).

Figure 1. Immune scores and stromal scores are associated with Kaplan-Meier survival in STAD. The OS and RFS curves between high and low groups based on immune scores (A, B). The OS and RFS curves between high and low groups based on immune scores (C, D).

The different distribution of stromal scores and immune scores in different clinicopathological features, including stage N (A, B) stage M (C, D) grade (E, F) and stage T (G, H). The different expression of CD274 between high and low groups based on immune scores and stromal scores (I, J). The different distribution of stromal scores and immune scores between different TP53 mutation and wildtype (K, L).

Figure 2. The different distribution of stromal scores and immune scores in different clinicopathological features, including stage N (A, B) stage M (C, D) grade (E, F) and stage T (G, H). The different expression of CD274 between high and low groups based on immune scores and stromal scores (I, J). The different distribution of stromal scores and immune scores between different TP53 mutation and wildtype (K, L).

Identification and analysis of DEGs and Hub genes

Next, the Agilent G microarray data of all 510 STAD cases obtained from the TCGA database was examined to find out the association of gene expression profiles. Heatmaps showed different gene expression profiles in the high immune-related scores vs low immune-related scores group (Figure 3A, 3B). Next, 223 DEGs were selected based on immune scores and stromal scores (fold change >1.0, P < 0.05), (Figure 4B). Also PCA and PPI network were constructed using the ClustVis (https://biit.cs.ut.ee/clustvis/), (Figure 3C, 3D and Figure 4A).

Heatmaps (A, B) and Principal Component Analysis (PCA) (C, D) showed different gene expression profiles in the immune scores of top half (high score) vs. bottom half (low score). P 1.5. Genes with higher expression are shown in red, lower expression are shown in blue.

Figure 3. Heatmaps (A, B) and Principal Component Analysis (PCA) (C, D) showed different gene expression profiles in the immune scores of top half (high score) vs. bottom half (low score). P <0.05, fold change >1.5. Genes with higher expression are shown in red, lower expression are shown in blue.

The KEGG pathway and GO enrichment analysis of DEGs were performed by DAVID (Figure 4C). The results of GO analysis revealed that the biological processes of DEGs showed significant enrichment in immune response, proteolysis, inflammatory response b and regulation of immune response. Molecular functions of DEGs demonstrated advancement in antigen binding, serine-type endopeptidase activity and transmembrane signaling receptor activity. The cell components that are enriched with DEGs included plasma membrane, integral component of membrane and extracellular region. KEGG pathway analysis revealed that DEGs were mainly enriched during cytokine-cytokine receptor interaction, chemokine signaling pathway, Staphylococcus aureus infection, and B cell receptor signaling pathway.

The PPI network of DEGs was constructed using Cytoscape. (A) Venn diagram was generated on the stromal scores and immunes scores of DEGs via Draw Venn Diagram (B). The biological process analysis of hub genes was constructed using BiNGO. The color depth of nodes refers to the corrected P-value of ontologies. The size of nodes refers to the numbers of genes that are involved in the ontologies. PC).

Figure 4. The PPI network of DEGs was constructed using Cytoscape. (A) Venn diagram was generated on the stromal scores and immunes scores of DEGs via Draw Venn Diagram (B). The biological process analysis of hub genes was constructed using BiNGO. The color depth of nodes refers to the corrected P-value of ontologies. The size of nodes refers to the numbers of genes that are involved in the ontologies. P<0.05 was considered statistically significant. The KEGG pathway and GO enrichment analysis of DEGs (C).

Eight hub genes were selected using BiNGO and were presented in Table 1. The selection criteria were as follows: MCODE scores >5, degree cut-off=2, node score cut-off=0.2, Max depth=100 and k-score=2. Next, Kaplan-Meier plots of the 8 hub genes were upregulated in the patients, decreasing the overall survival rate (Figure 5A). Interestingly, compared to normal tissues, GPNMB and TNFSF8 were significantly over-expressed in patients with STAD (P < 0.05) (Figure 5B).

Kaplan-Meier plots the 8 hub genes by R v3.5.1 from TCGA database. (A) The different expression levels of 8 hub genes expression between normal and tumor from TCGA (B).

Figure 5. Kaplan-Meier plots the 8 hub genes by R v3.5.1 from TCGA database. (A) The different expression levels of 8 hub genes expression between normal and tumor from TCGA (B).

Table 1. Functional roles of 8 hub genes with degree ≥10.

Gene symbleFull nameFunction
CCDC80Coiled-Coil Domain Containing 80Diseases associated with CCDC80 include Herpes Zoster Oticus and Peripheral Nervous System Disease. Gene Ontology (GO) annotations related to this gene include heparin binding and glycosaminoglycan binding..
FCN1Ficolin 1Among its related pathways are Innate Immune System and Creation of C4 and C2 activators. Gene Ontology (GO) annotations related to this gene include calcium ion binding and antigen binding.
GPNMBGlycoprotein NmbAmong its related pathways are Adhesion and Signaling by GPCR. Gene Ontology (GO) annotations related to this gene include heparin binding and integrin binding. An important paralog of this gene is PMEL.
OLR1Oxidized Low Density Lipoprotein Receptor 1Among its related pathways are Innate Immune System and Cell surface interactions at the vascular wall. Gene Ontology (GO) annotations related to this gene include carbohydrate binding and low-density lipoprotein particle receptor activity.
PDCD1LG2Programmed Cell Death 1 Ligand 2Involved in the costimulatory signal, essential for T-cell proliferation and IFNG production in a PDCD1-independent manner. Interaction with PDCD1 inhibits T-cell proliferation by blocking cell cycle progression and cytokine production.
RESSF2Ras Association Domain Family Member 2RASSF2 (Ras Association Domain Family Member 2) is a Protein Coding gene. Among its related pathways are Hippo signaling pathway - multiple species. An important paralog of this gene is RASSF4.
TNFSF8TNF Superfamily Member 87Diseases associated with TNFSF8 include Anaplastic Large Cell Lymphoma and Lymphoma, Hodgkin, Classic. Among its related pathways are ERK Signaling and Akt Signaling. Gene Ontology (GO) annotations related to this gene include signaling receptor binding and tumor necrosis factor receptor binding.
C2Synaptopodin 2 LikeActin-associated protein that may play a role in modulating actin-based shape.

Prognosis of GPNMB in STAD patients

A total of 317 STAD patients were obtained from the TCGA. The demographic and clinicopathological characteristics of patients are listed in Table 2. The results of multivariate analysis of the 8 genes are listed in Figure 6A. Among these genes, GPNMB was considered as significantly risky gene with HR > 1 and P < 0.05. Multivariate analysis of GPNMB expression, related scores and clinicopathological characteristics revealed that its expression, grade, stage and stromal scores were statistically significant factors for the progression of STAD (Figure 6B). In addition, all the input factors should be incorporated into the nomogram later as adjustment items. The prognostic nomogram by integrating all significant independent factors from the multivariate analysis for predicting the overall survival rate in the training cohort was shown in Figure 6C. The area under the curve of ROC in prognostic nomogram for overall survival prediction was 0.86 (95%CI, 0.80 to 0.92) (Figure 6D). Interestingly, compared to low scores group, GPNMB showed significant overexpression in high immune-related score STAD group (Figure 6E). Also expression of GPNMB showed positive correlation with stage and immune infiltration including CD8+T, macrophages, neutrophils and dendritic cells (Figure 6F, 6G).

Multivariate Cox regression analyses among the hub genes. (A) and related clinical characteristics respectively (B). The prognostic nomogram that integrated all significant independent factors from the multivariate analysis for OS in the training cohort (C). The RCO curve area evaluating the prognostic nomogram for OS prediction was 0.860 (D). The relate GPNMB expression between high and low immune scores (E). The different expression of GPNMB in STAD with different stage features (F). The correlation of GPNMB expression with immune infiltration level in STAD via TIMER (G).

Figure 6. Multivariate Cox regression analyses among the hub genes. (A) and related clinical characteristics respectively (B). The prognostic nomogram that integrated all significant independent factors from the multivariate analysis for OS in the training cohort (C). The RCO curve area evaluating the prognostic nomogram for OS prediction was 0.860 (D). The relate GPNMB expression between high and low immune scores (E). The different expression of GPNMB in STAD with different stage features (F). The correlation of GPNMB expression with immune infiltration level in STAD via TIMER (G).

Table 2. Demographic and clinicopathological characteristics of patients with STAD in TCGA database.

Demographic or CharacteristicsAll subjects (N=317)(%)Training cohort (N=167)Validation cohort (N=150)
Age at diagnosis
<=4927(8.5)13(8.0)12(8.0)
50-5982(25.9)43 (26.0)39(25.8)
60-69109(34.4)58(34.5)51(34.4)
70-79105(33.1)56(33.5)49(33.0)
>=8022 (6.9)12(7.0)10(6.9)
Gender
Female120(37.9)63(37.8)57(38.8)
Male197(62.1)104(62.2)93(61.2)
Stage
I42(13.2)23(13.5)19(13.2)
II101(31.9)53(32.0)48(32.3)
III139(43.8)73(44.0)66(43.6)
IV35(11.0)18(10.5)17(9.9)
Grade
Grade I7(2.2)3(2.0)4(2.1)
Grade II108(34.1)57(34.0)51(33.0)
Grade III202(63.7)107(64.0)95(63.0)
T
T115(4.7)8(4.7)7(4.4)
T263(19.9)34(20.0)29(20.6)
T3152(47.9)82(49.0)85(47.8)
T487(27.4)47(28.0)40(27.5)
M
M0295(93.0)152(91.0)143(94.0)
M122(7.0)12(7.5)10(6.9)
N
N070(22.1)37(22.0)33(22.2)
N183(26.2)45(27.0)38(26.8)
N269(21.8)35(21.0)34(20.2)
N366(20.8)35(21.0)31(20.8)
GPNMB expression
High160(50.5)85(51.0)75(50.0)
Low157(49.5)82(49.0)75(50. 0)
Stromal scores
High158(49.8)84(50.0)74(50.0)
Low159(49.2)83(50.0)76(50.0)
Immune scores
High160(50.0)84(50.0)74(50.0)
Low157(50.0)83(50.0)76(50.0)

GPNMB promotes proliferation of STAD cell lines

Firstly, RT-qPCR assay was performed to investigate the expression of GPNMB in different STAD cell lines. The results showed that GPNMB expression was markedly upregulated in MGC803 and BSG823 cell lines (Figure 7A). The expression efficiency of pc-GPNMB was measured by RT-qPCR assay and western blotting in both MGC803 and BSG823 cell lines. Next, GPNMB was overexpressed by transfecting GPNMB plasmid into MGC803 and BSG823 cell lines, indicating that transfection of pc-GPNMB was effective and can be used for subsequent research (Figure 7B, 7C). The effects of GPNMB upregulation on cell proliferation and apoptosis were further examined in MGC803 and BSG823 cell lines. MTS analysis showed that the proliferation ability of MGC803 and BSG823 cells with GPNMB overexpression was significantly higher than that of the mock cells (Figure 7D). Colony formation analysis showed that GPNMB overexpression caused a significant increase in the number of colonies in MGC803 and BSG823 cell lines (Figure 7E).

The GPNMB expression in different STAD cells. (A) The expression efficiency of pc-GPNMB was measured by RT-Qpcr (B) and WB assay (C) in MGC803 and BSG823 cells. MTT assays show that GPNMB overexpression increased cell proliferation in MGC803 and BSG823 cells (D). Colony formation assays indicate significantly increased the number of colonies in pc-GPNMB group (MGC803 and BSG823 cells) cells compared to NC group (E). The plots showing tumor growth measurements of MGC803-pc-GPNMB, MGC803 (control), BSG823- pc-GPNMB and BSG823 (control) cells is shown, and the mean tumor weights in each group (MGC803-pc-GPNMB, MGC803, BSG823-pc-GPNMB and BSG823) on day 30 is shown (F). IHC staining showing that cell proliferation (Ki67-positive) positively correlates with GPNMB expression levels (G). All the above experiments were repeated six times respectively. *P

Figure 7. The GPNMB expression in different STAD cells. (A) The expression efficiency of pc-GPNMB was measured by RT-Qpcr (B) and WB assay (C) in MGC803 and BSG823 cells. MTT assays show that GPNMB overexpression increased cell proliferation in MGC803 and BSG823 cells (D). Colony formation assays indicate significantly increased the number of colonies in pc-GPNMB group (MGC803 and BSG823 cells) cells compared to NC group (E). The plots showing tumor growth measurements of MGC803-pc-GPNMB, MGC803 (control), BSG823- pc-GPNMB and BSG823 (control) cells is shown, and the mean tumor weights in each group (MGC803-pc-GPNMB, MGC803, BSG823-pc-GPNMB and BSG823) on day 30 is shown (F). IHC staining showing that cell proliferation (Ki67-positive) positively correlates with GPNMB expression levels (G). All the above experiments were repeated six times respectively. *P < 0.05 versus control.

GPNMB promotes STAD growth in vivo

To further study the biological significance of GPNMB in STAD patients, MGC803-pcGPNMB and BSG823-pcGPNMB cells were subcutaneously injected and their corresponding controls into nude mice and monitored their tumor growth. Compared to the control cells, MGC803-pcGPNMB and BSG823-pcGPNMB cells demonstrated a significantly larger tumor size. In addition, the high GPNMB levels in STAD cells showed association with higher tumor growth rates and larger tumor weight (Figure 7F). Immunohistochemical experiments revealed that the Ki67 positive rate in high-GPNMB group was significantly higher than that of the control group (Figure 7G). Taken together, these data suggest that GPNMB promotes the growth of STAD tumors in vivo.

mir-30b-3p targets GPNMB

Finally, two sets of miRNAs related to GPNMB were selected from miRanda and miRDB databases and obtained miR-30p-3b as a possible target candidate (Figure 8A). The expression of miR-30p-3b between tumors and normal tissues was compared via ENCORI. The results showed that the expression of miR-30p-3b in tumors was lower than that of normal tissues (P < 0.05) (Figure 8B). miRNA-target coexpression between miR-30p-3b and GPNMB via ENCORI showed that miR-30p-3b expression was negatively correlated with GPNMB expression (Figure 8C). More importantly, miR-30p-3b might be combined with 3’ UTR of GPNMB (Figure 8D). Therefore, miR-30p-3b might be speculated to act as a potential key upstream negative regulator of GPNMB and might be related to cancer treatment. The expression levels of GPNMB between miR-30p-3b and control groups in MGC803 and MSG823 cell lines was evaluated by qRT-PCR, and the results showed that GPNMB was down-regulated in both cells (Figure 8E, 8F). To determine the role of miR-30p-3b in regulating GPNMB expression, MGC803 cells or mutant miR-30p-3b mimics were transfected to luciferase-labeled NEK2-3'UTR, and then were analyzed by luciferase reporter gene analysis. The results showed that miR-30p-3b mimic has significantly reduced the luciferase activity of GPNMB-3'UTR, while the mutant miR-30p-3b mimic did not inhibit the luciferase activity of GPNMB -3'UTR (Figure 8G). Furthermore, the expression levels of GPNMB in MGC803 cells transfected with miR-30p-3b or negative control (NC) were analyzed by qRT-PCR and western blotting, and the results showed that the expression of GPNMB with miR-30p-3b was significantly lower than that of negative control group (Figure 8H). Therefore, these results indicated that GPNMB was upregulated by miR-30p-3b and might have an impact on the prognosis of STAD patients.

GPNMB is a target for miR-30b-3p in STAD cells. Venn diagrams showing the number of potential miRNAs targeting the 3’UTR of GPNMB, as predicted by two databases, miRanda and miRDB. (A) The gene expression of miR-30b-3p between cancer and normal samples via ENCORI. (B) The coexpression between miR-30b-3p and GPNMB via ENCORI showed that miR-30b-3p expression is negatively correlated with GPNMB expression. (C) Sequences of miR-30b-3p and their potential binding sites in the 3’UTR of GPNMB is shown. (D) Quantitative real time PCR analyzing miR-30b-3p expression relative to MGC803 as internal control is shown. (E) Comparison of GPNMB expression in STAD cells transfected with miR-30b-3p mimic or negative control (NC) based on qRT-PCR (F), and western blotting (H), and the loading control for western blotting was GADPH. Analysis of luciferase activity from reporters containing the 3’UTR end of GPNMB in cells transfected with the miR-30b-3p mimic, miR-30b-3p mutation mimic (miR-30b-3p-mut) and negative control (NC) is shown. (G) All the above experiments were repeated six times respectively. *P

Figure 8. GPNMB is a target for miR-30b-3p in STAD cells. Venn diagrams showing the number of potential miRNAs targeting the 3’UTR of GPNMB, as predicted by two databases, miRanda and miRDB. (A) The gene expression of miR-30b-3p between cancer and normal samples via ENCORI. (B) The coexpression between miR-30b-3p and GPNMB via ENCORI showed that miR-30b-3p expression is negatively correlated with GPNMB expression. (C) Sequences of miR-30b-3p and their potential binding sites in the 3’UTR of GPNMB is shown. (D) Quantitative real time PCR analyzing miR-30b-3p expression relative to MGC803 as internal control is shown. (E) Comparison of GPNMB expression in STAD cells transfected with miR-30b-3p mimic or negative control (NC) based on qRT-PCR (F), and western blotting (H), and the loading control for western blotting was GADPH. Analysis of luciferase activity from reporters containing the 3’UTR end of GPNMB in cells transfected with the miR-30b-3p mimic, miR-30b-3p mutation mimic (miR-30b-3p-mut) and negative control (NC) is shown. (G) All the above experiments were repeated six times respectively. *P < 0.05 versus control.

Discussion

Increasing evidence has shown that the differences in the effectiveness of tumor immunotherapy can be attributed to inhomogeneities of tumor microenvironment [9, 10]. The tumor microenvironment is composed of tumor cells, tumor infiltrating immune cells and extracellular matrix (ECM), which together promote tumor growth, transformation and invasion, protect tumor cells from immune surveillance by the host, and cause tumors to develop resistance and provide the basis for tumor dormancy and metastasis [11]. Recent studies have revealed that immune microenvironment directly or indirectly affects tumorigenesis and development. The mechanism involves promotion of tumor angiogenesis, changing of tumor biological characteristics, screening of tumor cells that are adapted to the microenvironment for survival, and promotion of tumor progression by establishing a suitable tumor microenvironment [12]. In this study, we demonstrated that patients with low stromal scores had slightly increased overall survival and disease-free survival rates when compared to that of high stromal scores, while the overall survival rate of low immune scores was much higher than that of high immune scores. The stromal scores and immune scores were positively correlated with stages T and M and grade, indicating that the related immune scores might contribute to STAD progression (Figure 1B, P < 0.005).

Sun et al. have measured the expression levels of PD-L1 in cancer tissues and normal tissues of 102 patients with gastric cancer, and the results showed that the positive expression rate of PD-L1 in gastric cancer tissues was 40%-63%, while hardly detected in normal gastric tissues [13]. Zhang et al. have conducted meta-analysis of 10 clinical studies including 1901 patients with gastric cancer, and found that the positive expression rate of PD-L1 in gastric cancer was significantly higher than that in normal gastric tissues, yielding better immunotherapy effects [14]. In this study, the expression of PD-L1 with different immune scores and stromal scores revealed that the expression of PD-L1 in high score group was higher than that of low score group. Also the high expression of PD-L1 might induce apoptosis of anti-tumor T lymphocytes, and achieve immune escape of gastric cancer cells, promoting the occurrence and development of gastric cancer. Many studies have shown that TP53 mutations play a negative role in antitumor immunity [15, 16]. Nilay revealed that deletion of TP53 in gastric cells conferred a selective advantage and promoted the development of dysplasia in the setting of dietary carcinogenesis [17]. TP53 mutations frequently occur in STAD patients and are associated with unfavorable clinical outcomes. To explain these findings, the immune scores and stromal scores of gastric cancer cohorts between TP53 mutant state from TCGA project were analyzed. The results showed that TP53 mutant cases had lower immune and stromal scores when compared to those with wild-type TP53. TP53 mutation itself results in decreasing the immune activities in STAD and other cancer types, and immune cell infiltration and immune activities showed a positive association with survival prognosis in STAD [18].

Next, 223 DEGs were selected based on the immune scores and stromal scores. The KEGG pathway and GO analysis showed significant enrichment of DEGs during immune response and its regulation. The 8 hub genes were selected using BiNGO. Kaplan-Meier plots of the 8 hub genes were upregulated in STAD patients, decreasing their overall survival. Interestingly, compared to normal tissues, GPNMB was significantly over-expressed in STAD (P < 0.05). GPNMB was discovered in 1995 by Weterman et al. [19], which is a type I transmembrane glycoprotein, and forms a new signal transduction pathway with Melanocyte Inducing Transcription Factor (MITF), and might promote the development of human cancer. Some reports have suggested that GPNMB might be involved in the differentiation of tissue cells and metastasis of tumor cells [20, 21], and is associated with the occurrence and invasion of melanoma cells, glioma cells, breast [2224], colorectal cancer [25], and prostate cancer [26] cells. Tomihari confirmed that GPNMB in patients with melanoma has the ability to downregulate the activation of melanoma-reactive T cells, thereby allowing melanoma to evade immunologic recognition and destruction [27]. Identification of biomarkers that distinguish responders and nonresponders might improve the management of patients with cancer. GPNMB checkpoint differs from that of the PD1 signaling pathway in expression and inhibitory mechanisms, and GPNMB expression might regulate immune checkpoint inhibitors responsiveness [28]. These studies were consistent with the findings of our study with regard to the association of high expression of GPNMB in STAD with worse prognosis and advanced progression. Interestingly, we found that the expression of GPNMB was positively correlated with immune infiltration cells including CD8+T, macrophages, neutrophils and dendritic cells. In addition, multivariate analysis revealed that age, grade, stage, stromal scores and GPNMB were considered as statistically significant factors for STAD progression. A prognostic nomogram has integrated all significant independent factors obtained from the multivariate analysis for overall survival. Wang et al. have developed a significant prognostic nomogram for predicting the respective overall survival rates of STAD, and the C-Index for overall survival prediction was 0.707 [29]. However, our nomogram has 8 readily available pathological variables and achieved a C-index of 0.860, which was superior to that of other nomograms and larger than that of the traditional seventh AJCC staging system (0.860vs 0.661). This indicated that high expression of GPNMB acts as an independent risk factor and improved the current prognostic model. Furthermore, cell function experiments and xenograft tumors in vivo were performed to further verify the role of GPNMB in tumor progression, and the data analysis results were presented above.

MicroRNAs (miRNAs) are a wide range of regulatory non-coding RNAs (18-25 nt) in mammals and can either activate or inhibit the expression of target genes through post-transcriptional regulation. At present, the abnormally expressed miRNAs are always used as key factors in disease prediction, prevention and treatment [3032]. Based on the miRanda and miRDB databases, miR-30b-3p might upregulate and bind to 3’UTR of GPNMB, which was confirmed by fluorescent enzyme reporter gene experiment. MiR-30b-3p is either upregulated or downregulated in several types of cancers. Previous studies have revealed that miR-30b-3p is downregulated in primary prostate cancer (PCa) and metastatic castration resistant PCa and can directly inhibit androgen receptor and PCa cell proliferation [33]. The expression of miR-30b-3p is markedly decreased in hepatocellular cancer tissues and cells, showing positive correlation with higher overall survival rate [34]. The expression of miR-30p-3b in tumors was higher than that of normal tissues (P < 0.05). More importantly, miR-30p-3b might be combined with 3 ’UTR of GPNMB, and miR-30p-3b expression was negatively correlated with GPNMB expression, which was confirmed by fluorescent enzyme reporter gene experiment. Therefore, we speculated that hsa-miR-346 acts as a potential key upstream negative regulator of GPNMB and might assist in treating cancer.

In summary, the prognostic value of stromal scores and immune scores in STAD was confirmed, and the hub genes were selected and analyzed from DEGs. Importantly, GPNMB was frequently overexpressed in STAD and was associated with aggressive STAD, and high GPNMB expression also indicated poor prognosis and progression in patients with STAD. Next, a nomogram was developed to predict the overall survival of patients with STAD. Finally, miR-30b-3p was identified as an upstream regulator of GPFMB expression, which might be developed as a new treatment strategy in reducing the development of STAD.

Materials and Methods

Microarray data

The gene expression profiles and clinical information were downloaded from TCGA database, which is a landmark in cancer genomics program, and involves molecular characterization of over 20,000 primary cancer and matched normal samples spanning 33 cancer types. The stromal scores and immunes scores were calculated by an open source web tool ESTIMATE (https://bioinformatics.mdanderson.org/estimate/). Kaplan-Meier plots were then constructed to analyze the overall survival and disease-free of the stromal scores and immunes scores. Next, the stromal scores and immune scores in STAD with different clinicopathological features, including stage T/N/M, CD274 expression and TP53 mutation, were systematically compared.

Identification and analysis of DEGs

Heatmaps and principal component analysis (PCA) were constructed using R software. Differentially expressed genes (DEGs) were then identified using the limma package of the R statistical software (Fold change > 1.5 and adjusted P < 0.05 were set as the cutoff values). Venn diagram was generated based on the stromal scores and immune scores of DEGs via Draw Venn Diagram (http://bioinformatics.psb.ugent.be/webtools/Venn/). The protein-protein interaction (PPI) network was made using The Search Tool for the Retrieval of Interacting Genes (STRING). The KEGG pathway and GO enrichment analysis of DEGs were performed by DAVID (http://david.ncifcrf.gov). The top signaling pathways with false discovery rate (FDR) of <0.05, and -log FDR >1.5 were identified as significant cut-off values.

Hub genes analysis and survival curve

The hub genes were selected using BiNGO (Cytoscape). Kaplan-Meier plots were constructed to analyze the overall survival and disease-free survival rates of the 5 hub genes by R v3.5.1. The expression levels of hub genes between normal and turmor from the TCGA database were compared using the R software.

Prognostic nomogram

Multivariate Cox regression analyses were performed to determine the prognostic value of the hub genes and related clinical characteristics. The nomogram was built with potential risk factors (P <0.05) based on multivariate Cox analysis using the R software package. The predictive performance of the nomogram was measured by operating characteristic (ROC) curves. The relation between high and low scores of GPNMB expression was performed, and also the expression of GPNMB in STAD with different stage features was compared. The correlation of its expression with immune infiltration level in STAD was performed via TIMER (Tumor IMmune Estimation Resource), which is a web server used for comprehensive analysis of tumor-infiltrating immune cells (https://cistrome.shinyapps.io/timer/).

Cell culture and transfection

STAD cell lines were cultured in DMEM medium containing 10% fetal bovine serum under the conditions of 37C, 5% CO2, and 100% saturated humidity. Passive culture was performed when the cells confluent to 70% to 80%. The cells were cultured in 6-well plates and synchronized with serum-free medium for 24 hours. Further, the cells were divided into pc-NC groups and pc-GPNMB groups. Next, we transfected the blank plasmid and GPNMB-interfering plasmid following the instructions of the Lipofectamine 2000 transfection kit.

Quantitative real-time polymerase chain reaction and Western blotting

Qqt-PCR was used to detect the changes in the GPNMB expression. The aforementioned reverse transcription products were tested by Takara's SYBRPremixExTaqTM on an ABI7900 instrument using qPCR, and GAPDH expression was utilized as an internal reference. The following primers were used: GPNMB - TGCCAAGCGATTTCGTGATGT; GAPDH - ACCCACTCCTCCACCTTTGA, CTGTTGCTGTAGCCAAATTCGT. Cells from pc-NC group and pc-GPNMB group cells were collected, and total cell protein was extracted using cell lysate. Then the protein concentration was determined by a BCA kit. According to the results of protein concentration detection, 25uL of protein samples were subjected to Western blot detection.

Colony formation assay

Transfected MGC803 and BSG823 cells were seeded into 60 mm culture dishes and cultured for 14 days in the complete medium. Then, cells were fixed using methanol for 10 min and stained with 0.1% crystal violet (Sigma-Aldrich) for 15 min. Next, dishes were photographed and colonies with over 50 cells were counted.

In vivo xenograft tumor growth

Twelve nude mice were divided into an experimental and a control group (n = 6). Nude mice were kept in an environment of about 25 °C with a suitable humidity, and were provided with sufficient food and water during the administration. Experimental grouping: 5 × 106 cells transfected with pc-GPNMB in MGC803 cell (pc-GPNMB in BSG823 cell) were suspended in 0.1 mL of PBS and injected subcutaneously into nude mice. The control group was given the same amount of normal saline. The volume and weight of the xenograft tumors in the nude mice were then measured every ten days. Thirty days later, the nude mice were euthanized, and the transplanted tumors were removed. All animal experiments were performed as approved by the Committee of Huaihe hospital of Henan university.

Immunohistochemical analysis

Antibody system immunohistochemistry kits were purchased from Roche, Switzerland. According to the method of Dowset et al. [35], pale yellow, yellow or brown particles appeared in the nuclei of cells as positive expression. Readers blinded to the patient's pathological data, observed the expression of the whole film under a microscope. The expression levels of GPNMB protein were divided into control group and high expression groups.

miRNA database analysis

The potential miRNAs targeting GPNMB were downloaded from miRanda and miRDB databases. We performed survival analysis of miR-30d-3p in STAD downloaded from TCGA. The expression values of miRNAs from miRNA-seq data were scaled with log2(RPM+0.01). Co-expression analysis for miR-30d-3p and GPNMB were performed via ENCORI, which mainly focuses on miRNA-target interactions and is an open-source platform for studying the miRNA-ncRNA, miRNA-mRNA and RBP-mRNA interactions from CLIP-seq, degradome-seq and RNA-RNA interactome data (https://web.archive.org/web/20110222111721/http://starbase.sysu.edu.cn/).

Luciferase reporter assay

The GPNMB gene promoter was cloned by RT-q PCR and DNA fragments from the 3’-UTR of GPNMB inserted into the luciferase reporter vector pGL3. The STAD cells (10,000/well) were seeded in triplicate in 48 well plates. Then, the pGL3- GPNMB -3’UTR reporter plasmids (100 ng) plus 5 ng of pRL-TK renilla plasmid and increasing levels of negative control (NC), miR-30d-3p or mutant miR-30d-3p mimics were co-transfected into the MGC803 cells. Luciferase activity analysis was next performed to calculate the luciferase activity ratio of the reporter plasmid and the internal reference.

Statistical analysis

Statistical analysis was performed with SPSS software and R software (version 3.5.1). Factors were identified as significant at (PP < 0.1) in the univariate analysis. Comparisons between groups were performed using independent sample t-tests, and multiple group comparisons were performed using single-factor variance, and pairwise comparisons were performed using LSD Lt test. PP < 0.05 was considered statistically significant.

Abbreviations

STAD: Stomach adenocarcinoma; TCGA: The Cancer Genome Atlas; GPNMB: Glycoprotein non-metastatic melanoma protein B; PCA: principal component analysis; ROC: receiver operating characteristic; DEGs: Differentially Expressed Genes; Pca: Prostate cancer; FDR: False Discovery Rate; TIMER: Tumor IMmune Estimation Resource.

Author Contributions

KY, LW, JZ, CW and CX conceived and designed the study, which were proofed by CW, SL and CQ. SL and KY analyzed the data and wrote the manuscript. All authors read and approved the final manuscript.

Acknowledgements

The authors gratefully acknowledge contributions from the TCGA Network.

Conflicts of Interest

The authors declare no conflicts of interest.

Funding

The study was funded by the Science and Technology Research Project of Henan Province in 2018 (Grant no. 182102310372) and the Science and Technology Research Project of Henan Province in 2019 (Grant no. 192102310091).

References

  • 1. Bray F, Ferlay J, Soerjomataram I, Siegel RL, Torre LA, Jemal A. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. 2018; 68:394–424. https://doi.org/10.3322/caac.21492 [PubMed]
  • 2. Katai H, Ishikawa T, Akazawa K, Isobe Y, Miyashiro I, Oda I, Tsujitani S, Ono H, Tanabe S, Fukagawa T, Nunobe S, Kakeji Y, Nashimoto A, and Registration Committee of the Japanese Gastric Cancer Association. Five-year survival analysis of surgically resected gastric cancer cases in Japan: a retrospective analysis of more than 100,000 patients from the nationwide registry of the Japanese gastric cancer association (2001-2007). Gastric Cancer. 2018; 21:144–54. https://doi.org/10.1007/s10120-017-0716-7 [PubMed]
  • 3. Wagner AD, Syn NL, Moehler M, Grothe W, Yong WP, Tai BC, Ho J, Unverzagt S. Chemotherapy for advanced gastric cancer. Cochrane Database Syst Rev. 2017; 8:CD004064. https://doi.org/10.1002/14651858.CD004064.pub4 [PubMed]
  • 4. Sun Y. Tumor microenvironment and cancer therapy resistance. Cancer Lett. 2016; 380:205–15. https://doi.org/10.1016/j.canlet.2015.07.044 [PubMed]
  • 5. Whiteside TL, Demaria S, Rodriguez-Ruiz ME, Zarour HM, Melero I. Emerging opportunities and challenges in cancer immunotherapy. Clin Cancer Res. 2016; 22:1845–55. https://doi.org/10.1158/1078-0432.CCR-16-0049 [PubMed]
  • 6. Umansky V, Blattner C, Fleming V, Hu X, Gebhardt C, Altevogt P, Utikal J. Myeloid-derived suppressor cells and tumor escape from immune surveillance. Semin Immunopathol. 2017; 39:295–305. https://doi.org/10.1007/s00281-016-0597-6 [PubMed]
  • 7. Fernandes JV, Cobucci RN, Jatobá CA, Fernandes TA, de Azevedo JW, de Araújo JM. The role of the mediators of inflammation in cancer development. Pathol Oncol Res. 2015; 21:527–34. https://doi.org/10.1007/s12253-015-9913-z [PubMed]
  • 8. Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, Treviño V, Shen H, Laird PW, Levine DA, Carter SL, Getz G, Stemke-Hale K, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013; 4:2612. https://doi.org/10.1038/ncomms3612 [PubMed]
  • 9. Beatty GL, Gladney WL. Immune escape mechanisms as a guide for cancer immunotherapy. Clin Cancer Res. 2015; 21:687–92. https://doi.org/10.1158/1078-0432.CCR-14-1860 [PubMed]
  • 10. Sautès-Fridman C, Cherfils-Vicini J, Damotte D, Fisson S, Fridman WH, Cremer I, Dieu-Nosjean MC. Tumor microenvironment is multifaceted. Cancer Metastasis Rev. 2011; 30:13–25. https://doi.org/10.1007/s10555-011-9279-y [PubMed]
  • 11. Swartz MA, Iida N, Roberts EW, Sangaletti S, Wong MH, Yull FE, Coussens LM, DeClerck YA. Tumor microenvironment complexity: emerging roles in cancer therapy. Cancer Res. 2012; 72:2473–80. https://doi.org/10.1158/0008-5472.CAN-12-0122 [PubMed]
  • 12. Junttila MR, de Sauvage FJ. Influence of tumour micro-environment heterogeneity on therapeutic response. Nature. 2013; 501:346–54. https://doi.org/10.1038/nature12626 [PubMed]
  • 13. Sun J, Xu K, Wu C, Wang Y, Hu Y, Zhu Y, Chen Y, Shi Q, Yu G, Zhang X. PD-L1 expression analysis in gastric carcinoma tissue and blocking of tumor-associated PD-L1 signaling by two functional monoclonal antibodies. Tissue Antigens. 2007; 69:19–27. https://doi.org/10.1111/j.1399-0039.2006.00701.x [PubMed]
  • 14. Zhang M, Dong Y, Liu H, Wang Y, Zhao S, Xuan Q, Wang Y, Zhang Q. The clinicopathological and prognostic significance of PD-L1 expression in gastric cancer: a meta-analysis of 10 studies with 1,901 patients. Sci Rep. 2016; 6:37933. https://doi.org/10.1038/srep37933 [PubMed]
  • 15. Li L, Li M, Wang X. Cancer type-dependent correlations between TP53 mutations and antitumor immunity. DNA Repair (Amst). 2020; 88:102785. https://doi.org/10.1016/j.dnarep.2020.102785 [PubMed]
  • 16. Lyu H, Li M, Jiang Z, Liu Z, Wang X. Correlate the TP53 mutation and the HRAS mutation with immune signatures in head and neck squamous cell cancer. Comput Struct Biotechnol J. 2019; 17:1020–30. https://doi.org/10.1016/j.csbj.2019.07.009 [PubMed]
  • 17. Sethi NS, Kikuchi O, Duronio GN, Stachler MD, McFarland JM, Ferrer-Luna R, Zhang Y, Bao C, Bronson R, Patil D, Sanchez-Vega F, Liu JB, Sicinska E, et al. Early TP53 alterations engage environmental exposures to promote gastric premalignancy in an integrative mouse model. Nat Genet. 2020; 52:219–30. https://doi.org/10.1038/s41588-019-0574-9 [PubMed]
  • 18. Jiang Z, Liu Z, Li M, Chen C, Wang X. Immunogenomics analysis reveals that TP53 mutations inhibit tumor immunity in gastric cancer. Transl Oncol. 2018; 11:1171–87. https://doi.org/10.1016/j.tranon.2018.07.012 [PubMed]
  • 19. Weterman MA, Ajubi N, van Dinter IM, Degen WG, van Muijen GN, Ruitter DJ, Bloemers HP. Nmb, a novel gene, is expressed in low-metastatic human melanoma cell lines and xenografts. Int J Cancer. 1995; 60:73–81. https://doi.org/10.1002/ijc.2910600111 [PubMed]
  • 20. Singh M, Del Carpio-Cano F, Belcher JY, Crawford K, Frara N, Owen TA, Popoff SN, Safadi FF. Functional roles of osteoactivin in normal and disease processes. Crit Rev Eukaryot Gene Expr. 2010; 20:341–57. https://doi.org/10.1615/critreveukargeneexpr.v20.i4.50 [PubMed]
  • 21. Rose AA, Siegel PM. Osteoactivin/HGFIN: is it a tumor suppressor or mediator of metastasis in breast cancer? Breast Cancer Res. 2007; 9:403. https://doi.org/10.1186/bcr1791 [PubMed]
  • 22. Tse KF, Jeffers M, Pollack VA, McCabe DA, Shadish ML, Khramtsov NV, Hackett CS, Shenoy SG, Kuang B, Boldog FL, MacDougall JR, Rastelli L, Herrmann J, et al. CR011, a fully human monoclonal antibody-auristatin E conjugate, for the treatment of melanoma. Clin Cancer Res. 2006; 12:1373–82. https://doi.org/10.1158/1078-0432.CCR-05-2018 [PubMed]
  • 23. Kuan CT, Wakiya K, Dowell JM, Herndon JE 2nd, Reardon DA, Graner MW, Riggins GJ, Wikstrand CJ, Bigner DD. Glycoprotein nonmetastatic melanoma protein B, a potential molecular therapeutic target in patients with glioblastoma multiforme. Clin Cancer Res. 2006; 12:1970–82. https://doi.org/10.1158/1078-0432.CCR-05-2797 [PubMed]
  • 24. Rose AA, Pepin F, Russo C, Abou Khalil JE, Hallett M, Siegel PM. Osteoactivin promotes breast cancer metastasis to bone. Mol Cancer Res. 2007; 5:1001–14. https://doi.org/10.1158/1541-7786.MCR-07-0119 [PubMed]
  • 25. Mokarram P, Kumar K, Brim H, Naghibalhossaini F, Saberi-firoozi M, Nouraie M, Green R, Lee E, Smoot DT, Ashktorab H. Distinct high-profile methylated genes in colorectal cancer. PLoS One. 2009; 4:e7012. https://doi.org/10.1371/journal.pone.0007012 [PubMed]
  • 26. Tsui KH, Chang YL, Feng TH, Chang PL, Juang HH. Glycoprotein transmembrane nmb: an androgen-downregulated gene attenuates cell invasion and tumorigenesis in prostate carcinoma cells. Prostate. 2012; 72:1431–42. https://doi.org/10.1002/pros.22494 [PubMed]
  • 27. Tomihari M, Chung JS, Akiyoshi H, Cruz PD Jr, Ariizumi K. DC-HIL/glycoprotein nmb promotes growth of melanoma in mice by inhibiting the activation of tumor-reactive T cells. Cancer Res. 2010; 70:5778–87. https://doi.org/10.1158/0008-5472.CAN-09-2538 [PubMed]
  • 28. Chung JS, Ramani V, Kobayashi M, Fattah F, Popat V, Zhang S, Cruz PD Jr, Gerber DE, Ariizumi K. DC-HIL/gpnmb is a negative regulator of tumor response to immune checkpoint inhibitors. Clin Cancer Res. 2020; 26:1449–59. https://doi.org/10.1158/1078-0432.CCR-19-2360 [PubMed]
  • 29. Wang CY, Yang J, Zi H, Zheng ZL, Li BH, Wang Y, Ge Z, Jian GX, Lyu J, Li XD, Ren XQ. Nomogram for predicting the survival of gastric adenocarcinoma patients who receive surgery and chemotherapy. BMC Cancer. 2020; 20:10. https://doi.org/10.1186/s12885-019-6495-2 [PubMed]
  • 30. Blenkiron C, Miska EA. miRNAs in cancer: approaches, aetiology, diagnostics and therapy. Hum Mol Genet. 2007; 16:R106–13. https://doi.org/10.1093/hmg/ddm056 [PubMed]
  • 31. Farazi TA, Spitzer JI, Morozov P, Tuschl T. miRNAs in human cancer. J Pathol. 2011; 223:102–15. https://doi.org/10.1002/path.2806 [PubMed]
  • 32. Fuse M, Nohata N, Kojima S, Sakamoto S, Chiyomaru T, Kawakami K, Enokida H, Nakagawa M, Naya Y, Ichikawa T, Seki N. Restoration of miR-145 expression suppresses cell proliferation, migration and invasion in prostate cancer by targeting FSCN1. Int J Oncol. 2011; 38:1093–101. https://doi.org/10.3892/ijo.2011.919 [PubMed]
  • 33. Kumar B, Khaleghzadegan S, Mears B, Hatano K, Kudrolli TA, Chowdhury WH, Yeater DB, Ewing CM, Luo J, Isaacs WB, Marchionni L, Lupold SE. Identification of miR-30b-3p and miR-30d-5p as direct regulators of androgen receptor signaling in prostate cancer by complementary functional microRNA library screening. Oncotarget. 2016; 7:72593–607. https://doi.org/10.18632/oncotarget.12241 [PubMed]
  • 34. Gao D, Zhou Z, Huang H. miR-30b-3p inhibits proliferation and invasion of hepatocellular carcinoma cells via suppressing PI3K/Akt pathway. Front Genet. 2019; 10:1274. https://doi.org/10.3389/fgene.2019.01274 [PubMed]
  • 35. Dowsett M, Nielsen TO, A’Hern R, Bartlett J, Coombes RC, Cuzick J, Ellis M, Henry NL, Hugh JC, Lively T, McShane L, Paik S, Penault-Llorca F, et al, and International Ki-67 in Breast Cancer Working Group. Assessment of Ki67 in breast cancer: recommendations from the international Ki67 in breast cancer working group. J Natl Cancer Inst. 2011; 103:1656–64. https://doi.org/10.1093/jnci/djr393 [PubMed]