Comparing sequencing assays and human-machine analyses in actionable genomics for glioblastoma.

OBJECTIVE: To analyze a glioblastoma tumor specimen with 3 different platforms and compare potentially actionable calls from each.
METHODS: Tumor DNA was analyzed by a commercial targeted panel. In addition, tumor-normal DNA was analyzed by whole-genome sequencing (WGS) and tumor RNA was analyzed by RNA sequencing (RNA-seq). The WGS and RNA-seq data were analyzed by a team of bioinformaticians and cancer oncologists, and separately by IBM Watson Genomic Analytics (WGA), an automated system for prioritizing somatic variants and identifying drugs.
RESULTS: More variants were identified by WGS/RNA analysis than by targeted panels. WGA completed a comparable analysis in a fraction of the time required by the human analysts.
CONCLUSIONS: The development of an effective human-machine interface in the analysis of deep cancer genomic datasets may provide potentially clinically actionable calls for individual patients in a more timely and efficient manner than currently possible. CLINICALTRIALSGOV IDENTIFIER: NCT02725684.

The clinical application of next-generation sequencing technology to cancer diagnosis and treatment is in its early stages.[1-3] An initial implementation of this technology has been in targeted panels, where subsets of cancer-relevant and/or highly actionable genes are scrutinized for potentially actionable mutations. This approach has been widely adopted, offering high redundancy of sequence coverage for the small number of sites of known clinical utility at relatively low cost.

However, recent studies have shown that many more potentially clinically actionable mutations exist both in known cancer genes and in other genes not yet identified as cancer drivers.[4,5] Improvements in the efficiency of next-generation sequencing make it possible to consider whole-genome sequencing (WGS) as well as other omic assays such as RNA sequencing (RNA-seq) as clinical assays, but uncertainties remain about how much additional useful information is available from these assays.

Aside from cost, a challenge of WGS or whole-transcriptome data is the expertise and time required to interpret the full spectrum of somatic mutations. To address this challenge, Watson for Genomics (Watson Genomic Analytics [WGA]), a cancer analytic tool, uses standard variant call files (VCFs), copy number variant (CNV), and differential gene expression data to return a list of recommended cancer drugs. Here, we present the results of a targeted cancer panel along with WGS and RNA-seq in a patient with glioblastoma (GBM). We also compare results of expert interpretation of the tumor genome by bioinformaticians and oncologists at New York Genome Center (NYGC) and at collaborating institutions with those generated by WGA.

METHODS

Standard protocol approvals, registrations, and patient consents.

This study was approved by multiple institutional review boards (IRBs), including Rockefeller University IRB and Biomedical Research Alliance of New York IRB. The study was registered in ClinicalTrials.gov (NCT02725684). Informed written consent was obtained from the participant.

Participant.

This report describes the first participant in a multi-institutional study. NYGC-GBM-01 was a 76-year-old man with GBM. DNA and RNA were extracted from snap-frozen tissue. DNA from blood was obtained for comparison. The samples were analyzed by WGS and RNA-seq.

Single nucleotide variants and INDELs.

Whole-genome libraries were prepared using the Illumina TruSeq Nano DNA Sample Prep Kit and were sequenced on Illumina HiSeq X instruments (Illumina, San Diego, CA). Paired-end 2 × 150 bp reads were aligned to the GRCh37 human reference (BWA aln v.0.7.8)[6] and processed using a pipeline that includes marking of duplicate reads using Picard tools and realignment around INDELs and base recalibration using Genome Analysis Toolkit (GATK) version 2.7.4.[7] muTect v1.1.4,[8] LoFreq v2.0.0,[9] Strelka v1.0.13,[10] Pindel,[11] and Scalpel[12] were used to return the union of variant calls. Variants were filtered out if they were at >1% frequency in the 1000 Genomes or ExAC data sets, had more than 2 alleles to remove artifacts, raw frequency in the tumor was lower than that in the normal, or matched a custom “blacklist” of known systematic errors generated by comparing normal germline replicates. Remaining single nucleotide variants (SNVs) and INDELs were annotated via snpEff,[13] snpSift,[13] and GATK VariantAnnotator using annotation from ENSEMBL,[14] COSMIC,[15] Gene Ontology,[16] and 1000 Genomes.[17]

Structural variation.

Structural variants (SVs), such as CNVs and complex genomic rearrangements, were detected by NBIC-seq,[18] Delly,[19] CREST,[20] and BreakDancer.[21] We prioritized SVs in the intersection of callers and those with additional split-read evidence via SplazerS.[22] SVs with split-read support in the matched normal or annotated as known germline variants (1000 Genomes call set, Database of Genomic Variants) were removed as likely germline variants. The predicted somatic SVs were annotated with gene overlap (RefSeq, Cancer Gene Census) including prediction of potential effect on resulting proteins.

Tumor purity and ploidy.

Tumor purity was calculated from WGS data using Titan.[23] In addition, purity and ploidy were calculated from the Illumina OMNI 2.5M Array using ASCAT.[24]

RNA sequencing.

We used the Illumina TruSeq stranded messenger RNA protocol and sequenced 100 million reads. Reads were aligned using STAR[25] and Gencode genes were quantified using featureCounts.[26] Ninety-five percent of reads mapped the reference genome. We normalized the counts with DESeq2 and adjusted the quantification to account for GC bias[27] and batch effects[28] between The Cancer Genome Atlas (TCGA) GBM RNA-seq and our sample. The normalized expression data are used to identify GBM subtypes.[29]

Therapeutic targets and drug recommendations.

The NYGC uses the custom clinical Tier classification system for SNVs. Tier 1 variants are clinically important variants in the cancer type being studied (e.g., epidermal growth factor receptor [EGFR] T790M is known to be clinically important in lung cancer[30]). The same variant observed in a cancer unknown to manifest this variant is classified as Tier 2 (e.g., the clinical importance of EGFR T790M is unknown in GBM). Tier 3 variants are in targetable genes; however, the specific variant is not known to be targetable (e.g., an unknown mutation in EGFR). Tier 4 variants are in genes cataloged by COSMIC cancer census and not included in Tiers 1–3.[15] All other variants are in Tier 5 and considered variants of uncertain significance (VUS). Variants in Tiers 1–4 are considered potentially targetable. Variants were matched to potential treatments by identifying the most aberrant genes from a combination of SNV, INDEL, SV, and RNA-seq data and by searching the NYGC drug-to-gene database. Prioritization of potential treatments was based on further manual assessment including criteria such as strength of data supporting variants detected, FDA approval of drug in GBM or in another cancer type, current GBM trial for a drug, and successful use of the drug to target the variant identified to treat GBM or other cancer types.

Watson Genomic Analytics.

WGA, an IBM research proof-of-concept environment of Watson for Genomics,[31] is a cognitive system built on several different predictive models to analyze up to whole-genome scale molecular data. VCFs, CNV, and gene expression data are input to WGA. The VCF file provided to WGA contains the union from 3 calling algorithms each for SNVs and INDELs specified in the Methods section. CNV data are inputted as copy number log2 (T/N) ratio values per gene. Modified Z-scores of RNA-seq normalized expression data per gene are used as proxy for differential gene expression. Modified z-score per gene is calculated by subtracting the median transcripts per million (TPM) value (over the TCGA GBM cohort) to this sample's TPM and dividing by the TCGA SD. With this input, WGA leverages a comprehensive database of structured (20+ sources include DrugBank, NCI, COSMIC, ClinVar, and 1000 Genomes) and unstructured (evidence extracted from literature using Natural Language Processing [NLP]) biological and medical data. To date, WGA processed abstracts from PubMed and where possible, began analyzing full-text articles. In addition, the NLP engine is being trained to understand the approximate 5,600 clinical trials at ClinicalTrials.gov. It is from the unstructured sources that WGA maintains a current repository of drug-disease associations and biomarkers for prognosis and therapeutics, as well as matching patients to relevant clinical trials based on molecular criteria. WGA identifies gene alterations most likely to be important in cancer and then identifies relevant treatments that directly or indirectly target the variant. WGA also identifies VUS, resistive or sensitizing markers for the drug of interest, and relevant clinical trials.

RESULTS

Case report.

NYGC-GBM-01 was a 76-year-old man who presented with headache and difficulty with ambulation. CT of the brain revealed a mass in the left parietal region. He underwent initial resection for which pathology revealed a GBM, negative for the following: EGFR amplification by in situ hybridization fluorescence, EGFRvIII RNA expression, IDH1 R132H by immunostaining, 1p36/19q13 deletion by fluorescent in situ hybridization analysis, and MGMT methylation. Postsurgically, he had right-sided hemineglect and right/left confusion. He became somnolent and required a re-resection and ventriculoperitoneal shunt placement. Two months after initial resection, he completed radiation therapy with 40 Gy over 3 weeks with concurrent temozolomide 75 mg/m2 daily. He then completed 3 cycles of adjuvant temozolomide at 5 months after initial resection, after which progression was seen on MRI. He considered multiple options, but experienced functional decline, and was no longer trial eligible. Instead, he started on the first dose of bevacizumab and CCNU, 7 months after initial resection. He further declined and died 1 month later.

A sample from the initial resection was examined with a FoundationOne test, and a snap-frozen sample was received for sequencing in this study. DNA and RNA extraction, sequencing, and analysis required 7 weeks at which time a tumor board meeting was convened, including the treating oncologist, a neuro-oncologist, and bioinformaticians, to discuss the results of the analysis completed by the NYGC. A clinical report of the findings was subsequently issued to the oncologist. This tumor board meeting occurred after the completion of the first cycle of adjuvant temozolomide. The oncologist planned on referring the patient to clinical trials identified by the NYGC, but at the time of progression, he was no longer trial eligible due to functional decline.

Tumor analysis.

The metrics for the WGS are shown in table 1. The sample had an estimated tumor purity of 47%–52% and ploidy of 1.99. Table 2 describes the types and number of variants identified. Specifically, variant-calling pipeline analysis identified 8,449 total somatic mutations (with 150 falling in exonic, protein-coding regions) and a complex landscape of amplifications and deletions. Mutational signatures are an important molecular characterization of the tumor and assessing applicability of immunotherapy,[32] for example, RNA-seq identified the sample as the mesenchymal subtype of GBM.[28] WGS supported this with evidence of an NF1 mutation and a CDKN2A loss with a gain of Chr 7 and a loss of Chr 10.

Table 1

Sample whole-genome sequencing metrics

Table 2

Number of somatic variants identified

The NYGC identified 6 actionable SNVs, of which 2 were Tier 3 variants (MET R755fs and FGFR3 L49V) and 4 were Tier 4 variants (in STAG2, PIK3R1, NF1, and ERG, described in table e-1 at Neurology.org/ng). In addition, 5 CNVs were identified, of which 2 were in genes that had SNVs (table 3). CREST and Pindel identified a 299-bp intragenic deletion at the intron-exon junction of exon 11 in MET, as well as an amplification of MET (log2 CNV 3.64-fold tumor vs normal amplification). RNA-seq confirmed overexpression of MET (z-score 2.23) and an in-frame exon-skipping event (METex11), at an allele frequency of approximately 50% (figure 1). This observation is molecularly analogous to the skipping of exon 14 identified in lung adenocarcinoma and other cancer types.[33] Although METex11 is located in the extracellular domain, we hypothesized that this mutation could lead to an overactivation of MET and could be targeted by a tyrosine kinase inhibitor. By analogy with studies of MET(Δ7–8),[34] both would lead to the lack of transmembrane localization. We also noted that mislocalization of MET(Δ7–8) renders the variant not targetable using antibodies. However, MET-specific tyrosine kinase inhibitors could efficiently deactivate the kinase.

Table 3

List of variants identified as actionable by 3 different platforms

Figure 1

Sashimi plot representing the MET exon–skipping event

Red lines indicate exon coverage and exon junctions. Numbers in red indicate the number of reads supporting these junctions (for instance, 1,181 reads are split between exons 10 and 12). Only junctions with more than 100 reads are represented here.

We also identified a codon insertion in PIK3R1 (p.R562_M563insIle/c.1686_1688dupTAT), which is a regulatory protein that interacts with and inhibits the functional catalytic protein, PIK3CA. Activating mutations of PIK3CA are known cancer drivers, as are loss-of-function mutations of PIK3R1. Functional studies have shown that PIK3R1 amino acid D560 is involved in hydrogen binding with PIK3CA and is an essential amino acid in regulating the activity of the catalytic subunit. The mutation identified here is in the same helical inhibitory (iSH2) domain of PIK3R1. The variant binds but fails to inhibit PIK3CA, leading to enhanced cell survival, Akt activation, anchorage-independent cell growth, and oncogenesis.[35,36] Furthermore, analysis of the crystal structure (figure 2) supports the conclusion that this mutation would inhibit the functional interaction between the 2 proteins, specifically through N345 of PIK3CA, resulting in PIK3CA activation.

Figure 2

Overlay of genomic mutations on crystal structure of PIK3R1

New York Genome Center identified GBM somatic mutations in iSH2 domain of PIK3R1 (red). Known D560Y and N564D (orange) mutants in the PI3K regulatory domain (PIK3R1) fail to inhibit PIK3CA and lead to enhanced cell survival, Akt activation, anchorage-independent cell growth, and oncogenesis.[34] The R562-M563 insertion mutation in NYGC-GBM1 potentially alters the functional interaction within this region, specifically through N345 (green) of PIK3CA. The potential inability of catalytic regulation due to PIK3R1 mutation can translate to therapeutic targeting with a PIK3CA inhibitor such as BKM120 and therefore the identification of trial NCT01870726 in combination with MET inhibitor INC280.

A recent cell line study identified a synergistic effect between MET exon 14–skipping variants and a PIK3CA E545K oncogenic variant, in which a combination of an MET inhibitor and a PIK3CA inhibitor showed better sensitivity than single therapy.[37] The PIK3CA E545K variant also activates PIK3CA. Taken together, these findings led us to suggest combinatorial INC280 (MET inhibitor) and BKM120 (PIK3CA inhibitor) therapy for potential clinical consideration, and this suggestion would have made the patient eligible for a clinical trial assessing efficacy of this combination (NCT01870726).

Watson for Genomic analysis.

Data for NYGC-GBM-01 were input into WGA, which produced a report summarizing actionable variants and a list of associated drugs, including some based on a pathway target analysis. WGA identified 6 actionable alterations, 14 associated drugs, 9 VUS (including FGFR3), copy number losses in Chr 9 (focal), 10 (armscale), and 11 (focal), and gain on Chr 7 (armscale). Both the NYGC and WGA identified 5 actionable alterations (in genes NF1, MET, CDKN2A, CDKN2B, and PIK3R1; table 3). WGA reported an NF1 SNV and annotated the variant as inactivating but did not deem copy number change to be sufficient for calling this or EGFR and PTEN. A 1-copy gain of EGFR was below WGA's threshold for classification as a targetable variant. Furthermore, this variant was shown to be negative for amplification by in situ hybridization fluorescence. However, the NYGC decided to list it as potentially targetable, given it is a known actionable variant in GBM. Similarly, a 1-copy PTEN loss is reported by the NYGC but not by WGA or by the FoundationOne. The NYGC reported this variant because of its clinical implications; it is associated with resistance to EGFR tyrosine kinase inhibition via AKT/mTOR pathway activation and is linked to cetuximab resistance, but can be targeted by mTOR inhibitors.[38] For PIK3R1, the NYGC identified BKM120 as a potential therapeutic option based on additional RNA-seq evidence of overexpression of PIK3CA. WGA identified PIK3R1 as a relevant variant via SNV data by WGS and used RNA-seq information in pathway and drug analysis to also recommend BKM120. The MET amplification and associated drugs are reported by both platforms; however, WGA had 2 drugs for MET amplification, whereas the NYGC prioritized 1 therapeutic option, INC280, based on GBM trial data availability. The NYGC reported 8 clinical trials associated with 5 genes. WGA found 10 clinical trials that may be relevant across 6 actionable alterations.

Comparison with a panel.

Table 3 also compares all variants and drugs identified by the NYGC and FoundationOne. NYGC analysis identified 8 unique variants not found by the FoundationOne, including an exon-skipping event. The NYGC identified drugs for 10 targets, while FoundationOne identified drugs for 4. Furthermore, 6 of the variants reported as of unknown significance occurring in the tumor by FoundationOne were germline variants. One variant (DNMT3A splice site 2083-1G>C) was called by FoundationOne, but the position and base change were different from a nearby variant identified by the NYGC.

DISCUSSION

The NYGC is undertaking a WGS research study in patients with GBM to investigate the efficiency and feasibility of WGS to inform therapeutic options. Here, the results of NYGC WGS and RNA-seq were compared with a clinical panel assay. Also, in collaboration with IBM, the NYGC examined the therapeutic options identified by WGA based on WGS and RNA-seq data. Genomic results from this patient clearly displayed the diversity of driver events typically seen in GBM. Of interest, we identified mutations in targetable genes that were not precise matches to known specific targetable variants, and which nonetheless suggested potential therapeutic options.

Multimodal analysis (WGS and RNA-seq) increased confidence in the identification of the MET mutation; analysis of the literature of prior MET exon–skipping events suggested the plausibility of considering a tyrosine kinase inhibitor that could target MET. Similarly, manual literature search of the PIK3CA E545K oncogenic variant led to the conclusion that this was likely an activating mutation. Moreover, manual database searches resulted in the suggestion of a combinatorial treatment with an MET inhibitor and a PIK3CA inhibitor, which made logical sense and also made the patient eligible for a clinical trial for this combination (NCT01870726).

None of these observations were evident from the panel. This suggests that pursuing a more extensive comparison of panel and deeper sequencing (e.g., WGS and RNA-seq) will be of interest. An added point, previously noted by others,[39] is that the sequencing of both germline and tumor DNA not only heightened our sensitivity for what variants might be tumor drivers but was able to rule out a number of germline variants called by the FoundationOne as not likely to be primary drivers of this patient's GBM.

Although we conducted WGS of this sample at roughly twice the cost of WES, the primary analysis was performed on the protein-coding region of the genome. There may be technical advantage to WGS even for assaying targeted regions. WES relies on hybridization capture of specific genes which introduces intrinsic bias for each gene as a function of GC/AT content, while WGS relies more simply on mechanical shearing of DNA prior to sequencing. Previous studies have found that for disorders caused by constitutional mutations,[40,41] WGS is more sensitive than WES for variant detection. To assess whether WGS could detect variants not identified by WES to justify the added cost, it would require a direct comparison of the assays on the same sample. We are undertaking a study to address this question.

This patient died approximately 8 months from the time of initial resection falling short of the median survival time for GBM. The oncologist recommended enrollment in a clinical trial targeting PIK3 and MET alterations on recurrence on adjuvant temozolomide. However, the patient's clinical decline eliminated his ability to participate in trials. This highlights one of the challenges of the clinical application of precision medicine technology. The identification of targets and potentially useful drugs in a timely manner is only the first step. Drug and drug trial access is crucial to determine the benefit of this approach in cancer management.

Another key observation was that the WGA analysis vastly accelerated the time to discovery of potentially actionable variants from the VCF files. As previously reported, we found that WGA was able to provide reports of potentially clinically actionable insights within 10 minutes, while human analysis of this patient's VCF file took an estimated 160 hours of person-time. This is critical if sequencing is to be brought out of the research arena and into the scaled, real-world clinical realm. This study is an important step forward promoting human-machine interface as a way to address a key bottleneck in cancer genomics.