Risk stratification of gastrointestinal stromal tumors by Nanostring gene expression profiling

The risk assessment classification schemes for gastrointestinal stromal tumors (GIST) include tumor site, size, mitotic count and variably tumor rupture. Heterogeneity in high-risk GIST poses limitations for current classification schemes. This study aims to demonstrate the clinical utility of risk stratification by gene expression profiling (GEP) using Nanostring technology. Fifty-six GIST cases were analyzed using a 231 gene expression panel. GEP results were correlated with clinical and pathological data. The prognostic performance was assessed in 34 patients with available survival data using ROC curves, Kaplan–Meier survival curves and compared with traditional risk assessment schemes. Volcano plot analysis identified seven genes with significantly higher expression (FDR < .0.05) in high-risk than in non-high-risk tumors, namely TYMS, CDC2, TOP2A, CCNA2, E2F1, PCNA, and BIRC5. Together, these transcripts exhibited significantly higher expression in high-risk tumors than in intermediate (P < 0.01), low (P < 0.001), and very low (P = 0.01) risk tumors. Receiver-operating characteristic curve analysis demonstrated area under the curve (AUC) to be 0.858 for the separation of high-risk and non-high-risk tumors. Kaplan–Meier survival analysis demonstrated improved risk stratification (log-rank test P < 0.001) compared to the current risk assessment classification (P = 0.231). In addition to current clinical and histology-based risk classification for patients with GIST, gene expression may offer complementary prognostic information.


Introduction
Gastrointestinal stromal tumors (GISTs) are the most common mesenchymal tumors of the gastrointestinal tract. GISTs show equal distribution across geographic regions and ethnic groups, and have slightly male (55%) or possibly no gender predisposition and annual global incidence of 10-15 cases per million (Miettinen et al. 2005;Soreide et al. 2016).
GISTs arise from interstitial cells of Cajal. New evidence suggests that PDGFRA-mutant GISTs are derived from telocytes rather than typical ICC cells (Ricci et al. 2018).
A number of current schemes for risk assessment in newly diagnosed primary GIST include three factors: mitotic index, tumor size, and tumor location Schmieder et al. 2016). A study using these factors showed that the median overall survival of all patients with diagnosed GIST was 11.7 years (Blanke et al. 2008). Some models for prediction of risk of recurrence, such as the National Institutes of Health (NIH) classification, consider tumor size and mitotic count alone (Fletcher et al. 2002). Studies using this criteria estimate the risk of recurrence in high-risk GISTs at 38.5% (Nakamura et al. 2005), 46% (Tryggvason et al. 2005), and 62.5% (Nilsson et al. 2005). Alternative methods incorporate other criteria, such as tumor site and tumor rupture Joensuu et al. 2014). A study using this method showed that as many as 76.9% patients with the highest risk showed recurrence during follow-up (Joensuu et al. 2014). These risk prediction models have been proposed as instruments for determining the potential benefits of adjuvant therapy (Huang et al. 2007). However, prognostic heterogeneity in high-and intermediate-risk patients, especially in the NIH scheme, is one of their limitations (Huang et al. 2007).
Gene expression profiling (GEP) is a powerful tool for assessing risk recurrence in several cancers (Gnant et al. 2014). Numerous GEP Studies in GIST have contributed to our knowledge of the pathobiology of GIST. Microarray analysis has shown distinctive profiles for gastric and intestinal GIST (Hara et al. 2015). AURKA was identified as a marker for high-risk GIST by bioinformatics analysis of publicly available GEP data and integration with clinicopathologic studies (Yen et al. 2012). CD133 was identified to be a marker for exon 11 mutation, gastric origin and poor prognosis by GEP (Arne et al. 2011). Hierarchical clustering of GEP results could be used to further subclassify GIST into two groups with distinctive clinical behavior and metastatic risk (Skubitz et al. 2016). Most of these studies used DNA microarray technique, a method that has not be easily translated into routine clinical laboratory workflow.
In this study, we hypothesized that gene expression signatures can be used to accurately classify patients with GIST into risk categories using Nanostring, a highly clinically reproducible GEP technique that is amenable to a routine clinical laboratory. The study furthers the goal of precision medicine by identifying and applying molecular determinants of GIST biological behavior to risk stratification.

Study cohort and clinico-pathological data
This study was approved by the Ethics Review Board at the University of Alberta (HREBA.CC-17-0375). The laboratory information system at the University of Alberta was searched for all cases with an original diagnosis of GIST between January 2010 and December 2016. In total, 56 GIST cases with residual FFPE material were identified and retrieved. In all cases, the initial diagnosis of GIST had been confirmed using CD117 and/or DOG1 immunohistochemistry. Pathological features, including tumor size, site, histological features, mitotic rate, risk assessment according to the modified NIH classification criteria, and patient tumor-node-metastasis (TNM) stage were extracted from the original pathology reports. All cases were reviewed and pathologic assessment confirmed by a soft tissue pathologist. Clinical follow-up data, which were available in the local cancer registry for 36 patients, were also retrieved. These clinical data included age at diagnosis, gender, postdiagnosis survival, and recurrence.

RNA isolation
One hematoxylin and eosin (H & E)-stained section and five unstained slides with 10-μm sections were obtained from each archival FFPE block. The H & E-stained sections were reviewed by a pathologist and tumor macrodissection was performed. RNA was isolated from the macrodissected sections using the RNeasy FFPE kit (Qiagen, Germantown, MD, USA). RNA concentration and purity were measured using a NanoDropTM 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA).

Gene expression analysis
A commercial gene expression panel, the NanoString GX Human Cancer panel (NanoString Technologies, Seattle, WA, USA), which includes 231 genes previously associated with human cancer, was used. These included transcripts associated with cell cycle, differentiation, proliferation, apoptosis, angiogenesis, and immune response. The panel was quantified using an nCounter FLEX system (NanoString Technologies, Seattle, WA, USA) as described previously (Mootha et al. 2003). Raw gene expression counts were quality controlled and normalized using the nSolver analysis software version 4.0 (NanoString Technologies, Seattle, WA, USA) and the positive and negative controls provided by the manufacturer. Data were normalized to six housekeeping genes included in the commercial panel. The housekeeping genes were CLTC, GAPDH, GUSB, HPRT1, PGK1 and TUBB.

Data analysis
Post-normalization statistical analysis and visualization were performed with R version 3.4.3 (R Foundation for Statistical Computing, Austria) and gene set enrichment analysis (GSEA) version 3.0 (Broad Institute, Inc. Massachusetts Institute of Technology, Cambridge, Massachusets, USA, and Reagents of the University of California, San Diego, California, USA) (Mootha et al. 2003;Subramanian et al. 2005). Differences between patient groups were assessed using the Mann-Whitney U test for continuous data and Fisher's exact test for categorical data. Normalized transcript counts were used for individual gene analysis. Geometric means of normalized transcript counts were used for aggregate gene set analysis. Exploratory analysis was performed using heat maps with unsupervised hierarchical clustering by Euclidean distance. Differential gene expression between high-risk and non-high-risk tumors was assessed using volcano plot analysis incorporating log2 fold change and Student's t test P values, with correction for multiple comparisons using false discovery rate (FDR) adjustment. Receiver-operating characteristic (ROC) curve analysis was used to assess individual gene and aggregate gene set performance. Cox proportional hazards regression was used to assess the prognostic significance of clinical, pathological, and molecular data. Kaplan-Meier curves were used to visualize differences in patient survival between patient groups. P values < 0.05 were considered statistically significant.
Biological pathway analysis was performed using GSEA v3.0. Gene Ontology (Ashburner et al. 2000) and Reactome (Fabregat et al. 2018) were used as the gene set databases. Parameters were set at 1000 phenotype permutations. FDR < 0.25 were used for pathway discovery. FDR < 0.05 was considered significant.

Immunohistochemistry (IHC) validation of GE results
PCNA and Ki67 IHC of tumor was performed and analyzed using the software QuPath 2.3 (Bankhead et al. 2017). The percentage of positive nuclei (PPN) of only malignant cells was computed using a trained cell classification algorithm applied to whole slide images (WSI). Receiver-operating characteristic (ROC) curves were constructed to assess performance for IHC risk prediction when correlated with pathological data. Kaplan-Meier survival curves were plotted using available patient survival data to compare prognostic performance with actual patient outcome.
Mutations in KIT and PDGFRA were analyzed for all samples during initial clinical workup by Sanger sequencing (Table 1). Forty-eight cases (86%) had detectable KIT mutations, whereas PDGFRA mutations were detected in only 5 cases (9%). Three cases (5%) had no detectable KIT or PDG-FRA mutations. These three cases were excluded from GEP modeling due to low number in this group. The detectable KIT mutations were located in exon 9 in 2 case (2%), exon 11 in 43 cases (77%), exon 13 in 1 case (2%), and exon 17 in 2 cases (4%). The PDGFRA mutations were detected in exon 18 in all 5 cases. The mutation types included singlenucleotide variants, insertion-deletions, duplications, and frame-shift changes.

Differential gene expression analysis
Exploratory analysis with hierarchical clustering identified two specific molecular signatures involving the top 50 differentially expressed genes. These signatures were capable of delineating the high-risk and non-high-risk groups as well as the low-risk and non-low-risk groups (Supplemental Figs. 1,2). Volcano plot analysis identified seven genes with statistically significant higher expression (FDR < 0.05) in the high-risk group than in the non-high-risk groups, namely TYMS, CDC2, TOP2A, CCNA2, E2F1, PCNA, and BIRC5 (Fig. 1). As an aggregate seven-gene set, these transcripts exhibited significantly higher expression in high-risk tumors than in intermediate (P = 0.002), low (P < 0.001), and very low (P = 0.007) risk tumors (Fig. 2). ROC curve analysis demonstrated an area under the curve (AUC) of 0.858 (95% CI 0.762-0.954) for the separation of high-risk vs. non-risk tumors using the seven-gene set expression (Fig. 3). This was superior to that observed for any of the individual genes (AUC, 0.770-0.840).

Survival analysis
To investigate whether the expression of the seven-gene set was a better predictor of outcome than the current clinicopathological risk assessment classification, Cox proportional hazards regression was performed for the 34/56 patients with available clinical follow-up data, among which 3 died during the study follow-up period. Univariate analysis did not demonstrate any significance in all the characteristics tested. This is likely due to our low sample size. We note, however, that the seven-gene set expression was the only variable very close to statistically significant association of increased risk of death (P = 0.05) ( Table 2). None of the other clinical or pathological features currently used for risk stratification, including tumor site, tumor size, and mitotic activity were statistically significant in this study (Table 2). We considered the results of the seven-gene set to be suggestive because it performed better than other characteristics that are known to be significant in other better powered studies. It is noteworthy that exon 11 mutations were present in all three patients who died although the hazard ratio for this variable could not be adequately analyzed owing to small sample size.
Kaplan-Meier survival analysis was performed to further assess the prognostic performance of the seven-gene set expression compared to the current risk assessment classification system. Interestingly, although gene selection for the seven-gene set was trained on the risk assessment labels of the existing classification system, statistically significant prediction of patient survival was only observed with the expression of the seven-gene set (log-rank test, P < 0.001) and not with the current risk assessment classification (P = 0.231) (Fig. 4).

Biological pathway analysis
To determine the specific biological pathways that determine risk stratification groups, GSEA was performed on the risk-stratified groups. The samples were classified into high-, intermediate-, and low/very low-risk groups. Pairwise comparison of risk groups was performed using GSEA. The pairs were high vs. non-high risk, intermediate vs. non-intermediate risk, and low/very low vs. non-low-risk. Each of the paired groups demonstrated differential up-and down-regulation of pathways. The high-risk group showed the most significantly divergent differential pathway enrichment from the non-high-risk groups (i.e., the intermediateand low/very low-risk groups). Four pathways from the Reactome gene sets were enriched in the high-risk group compared to non-high-risk groups, including cell cycle (FDR = 0.043), cell cycle regulation (FDR = 0.088), mitotic, G1, and S phases (FDR = 0.096), and interleukin signaling (FDR = 0.178). The seven-gene set previously identified using volcano plots of the differential expression analysis were also the most significant genes within the biological pathways identified (Supplemental Figs. 3, 4). These seven genes were frequently up-regulated in pathways enriched in the high-risk group and down-regulated in pathways in the low-risk group. Pathways with down-regulated genes in the high-risk group included developmental biology and axon guidance, although with FDR values of 0.42 and 0.37, respectively. One hundred and two pathways were enriched in the low-risk group compared to the non-low-risk group.
GIST gastrointestinal stromal tumors a Survival data were only available for 34/56 patients overall and 6/11 patients with high seven-gene set expression b Count greater than 71 The top six pathways from the Reactome gene sets were cell fate commitment, regulation of locomotion, regulation of GTPase activity, transmembrane receptor protein kinase activity, transmembrane receptor protein tyrosine kinase activity, and positive regulation of ERK1 and ERK2 cascade signaling (FDR < 0.25). The pathways identified in low-risk group can be generally categorized as pathways involved in cell differentiation and maturation and cell structure and function.
Leading edge analysis of these biological pathways was performed to identify the most significant genes within each pathway and the common genes among the pathways. Analysis showed a higher clustering of genes for high-risk groups than for low-risk groups. In the high-risk group, CDK1 was most highly expressed with PCNA, and was involved in all the top four biological pathways noted above (Supplemental Fig. 3). Other significantly expressed genes involved in four or more pathways included BRCA1/2, TOP2A, and E2F1. In the low-risk group, significantly expressed genes involved in the top four or more pathways included KDR, NTRK3, PDGFRB, EGFR, and ERBB4 (Supplemental Fig. 4).

IHC validation of seven-gene results
Since, the seven-gene sets were identified as playing a role in cell cycle signaling and progression, we chose two important gene markers of these process, namely Ki67 and PCNA for result validation by IHC. We performed image analysis of Ki67 and PCNA for risk groups and constructed ROC curves for Ki67 and PCNA, singly and in aggregate for risk prediction using the PPN score. Of the 34 cases with survival data, Ki67 and PCNA, alone and in combination, showed increased PPN in high-risk GISTS in comparison to non-high-risk GISTs (Figs. 5, 6). For Ki67, non-highrisk GISTs demonstrated a mean of 2.2% (range 0.7-4.5%) versus 14.0% (range 1.6-46.5%) for high-risk GISTs (P = 0.002), For PCNA, non-high-risk GISTs showed a mean of 13.2% (range 2.7-30.8%) versus 25.9% (range 3.6-83.5%) for high-risk GISTs (P = 0.041). Whereas, a combination of the two IHCs (product of PPNs) showed a mean of 32.2 (range 3.2-138.1) for non-high-risk GISTs and 570.0 (range 7.6-3881.5) for high-risk GISTs (P = 0.021). ROC curves were plotted for Ki67 and PCNA, singly and in combination to determine sensitivity, specificity, and to estimate the ideal cutoff for PPN. Ki67 alone was the most sensitive (sens 0.905, NPV 0.818). Whereas, the combined PCNA and Ki67 PPN was the most specific (spec 0.900, PPV 0.941). The optimal diagnostic cutoffs for PPN, as defined by Youden's J-statistic, were 2.782 for Ki67, 7.606 for PCNA, and 33.013 for Ki67 and PCNA combined. Performance of Ki67 and/ or PCNA PPN as potential risk prediction surrogates was assessed using Kaplan-Meier survival curves (Fig. 7). The probability of survival of high PPN cases was compared with low PPN cases. All three classes (Ki67 alone, PCNA alone, and combined Ki67 and PCNA) revealed poorer survival outcomes in specimens with high PPN in contrast to low PPN (Log-rank test, P < 0.0001). Of the high PPN pattern, combined Ki67 and PCNA demonstrated the worst overall survival outcomes. PCNA alone was the least predictive of poorer survival outcome.

Discussion
In this study, we aim to evaluate a potential role for gene expression analysis in risk stratification for GIST and to identify biological determinants for aggressive behavior. Using differentially expressed gene analysis, we identified seven genes, whose expressions signify poorer outcomes. We demonstrate their utility for stratification in high-risk GIST tumors.
Current risk stratification schemes in GIST utilize clinico-pathological criteria, including tumor size, site, mitotic count, and variably tumor rupture to predict clinical outcomes. The limitation of these methods is that a percentage of members in each risk group may show different biological behavior than the average member. Nanostring technology applies the principles of precision medicine to classify risk according to molecular features, which may predict behavior more accurately and can be used in a routine clinical setting.
We examined all GIST cases diagnosed at our institution from 2010 to 2016. Our study cohort showed an epidemiological distribution similar to those of other published cases in terms of median age, gender, site of occurrence, and KIT mutations. In this study, we analyzed various genes involved in cancer-related pathways, including cell cycle, cell differentiation, cell proliferation, apoptosis, angiogenesis, and immune response. Differential gene expression analysis identified molecular signatures for high-and lowrisk groups.
Seven genes, namely BIRC5, PCNA, E2F1, CCNA2, TOP2A, CDC2, and TYMS, demonstrated statistically significant differential expression (FDR < 0.05) between high-risk and non-high-risk tumors. These seven genes are involved in cell cycle control, apoptosis, and DNA replication. BIRC5 is a member of the inhibitor of apoptosis (IAP) gene family, which encodes negative regulatory proteins that prevent apoptotic cell death. PCNA is a cofactor of DNA polymerase delta, which assists in increasing the processivity of leading strand synthesis during DNA replication. E2F1 is a member of the E2F family of transcription factors and controls the cell cycle and the tumor suppressor protein activity. CCNA2 belongs to the cyclin family and is a regulator of the cell cycle. It mediates transition through the G1/S and G2/M checkpoints. TOP2A encodes a DNA topoisomerase and is involved in chromosome condensation, chromatid separation, and relaxation of DNA torsional stress. TYMS encodes the enzyme thymidylate synthase, which maintains the dTMP pool necessary for DNA replication and repair. Many of these pathways have been implicated by previous studies in GIST progression, including cell cycle control (Nakamura et al. 2005). Further analysis using GSEA confirmed a role for additional pathways such as cell cycle control and regulation, mitosis, and interleukin signaling. These pathways showed up-regulation of gene expression in high-risk GIST. A larger number of pathways were upregulated in low-risk and down-regulated in non-low-risk GIST. These pathways involve genes involved in cell fate commitment, cell differentiation and maturation, and cellular structure and function. This observation supports a biological mechanism of aggressive tumor behavior in GIST, in which a larger tumor fraction is in the mitotic phase of the cell cycle. This is consistent with the higher mitotic count seen in current pathologic assessment of high-risk GIST. The genes identified by our study were largely different from 25 genes previously identified as prognostic/predictive markers for soft tissue sarcomas, except for the single overlap of CDK1 (Takahashi et al. 2014).
We observe that the expression of the seven-gene set may be a better predictor of overall survival than traditional clinico-pathological features used for risk stratification, including age, gender, site of occurrence, tumor size, and mitotic rate especially for high-risk GIST. However, we noted that this analysis is limited by the small number of deaths and overall small number of samples tested in this series, especially for the low-risk group. For instance, all patients who died harbored an exon 11 mutation, but the mathematical requirement for a non-zero denominator in calculating the hazard ratio limited its assessment.
We validated our results with immunohistochemistry and image analysis, demonstrating that assessing markers of pathways involved by the seven-gene set was a useful predictor of disease risk especially in aggregate. While a comprehensive validation of all seven genes individually by IHC is ideal, for easy adaptation into routine histopathology use, we aimed for an IHC combination with two properties: simplicity of use for diagnostic translation;  Of note, all seven genes have important roles in various aspects of cell cycle activity. PCNA and Ki67 have been previously identified in various studies as markers of cell cycle activity and of increasing risk in GIST (Gumurdulu et al. 2007). On this basis, we selected these two genes to validate the performance of our seven-gene set.
Despite the limitations of this study, these results suggest that compared to the current clinical and histologybased risk assessment classification for patients with GIST, gene expression may offer significant but complementary prognostic performance. An arising question is whether heterogeneity in low-and intermediate-risk tumors can be defined by molecular markers. This interesting question could be addressed in future prospective studies with a larger cohort using our seven gene markers. Finally, Nanostring technology is simple, cost-effective and amenable to FFPE tissue, with established use in routine clinical GEP assays, which up till now has been hindered by inherent limitations of other techniques. Further validation studies are warranted to confirm the clinical utility of this molecular risk stratification approach.