Longitudinal analysis of treatment-induced genomic alterations in gliomas

Background Glioblastoma multiforme (GBM) constitutes nearly half of all malignant brain tumors and has a median survival of 15 months. The standard treatment for these lesions includes maximal resection, radiotherapy, and chemotherapy; however, individual tumors display immense variability in their response to these approaches. Genomic techniques such as whole-exome sequencing (WES) provide an opportunity to understand the molecular basis of this variability. Methods Here, we report WES-guided treatment of a patient with a primary GBM and two subsequent recurrences, demonstrating the dynamic nature of treatment-induced molecular changes and their implications for clinical decision-making. We also analyze the Yale-Glioma cohort, composed of 110 whole exome- or whole genome-sequenced tumor-normal pairs, to assess the frequency of genomic events found in the presented case. Results Our longitudinal analysis revealed how the genomic profile evolved under the pressure of therapy. Specifically targeted approaches eradicated treatment-sensitive clones while enriching for resistant ones, generated due to chromothripsis, which we show to be a frequent event in GBMs based on our extended analysis of 110 gliomas in the Yale-Glioma cohort. Despite chromothripsis and the later acquired mismatch-repair deficiency, genomics-guided personalized treatment extended survival to over 5 years. Interestingly, the case displayed a favorable response to immune checkpoint inhibition after acquiring mismatch repair deficiency. Conclusions Our study demonstrates the importance of longitudinal genomic profiling to adjust to the dynamic nature of treatment-induced molecular changes to improve the outcomes of precision therapies. Electronic supplementary material The online version of this article (doi:10.1186/s13073-017-0401-9) contains supplementary material, which is available to authorized users.

We have performed a deeper coverage of tumors as compared to matching blood samples (average target coverage was 194.3 and 121.3, respectively). The average percentage of reads with at least 20x coverage was 91.0% and 88.4% for tumor and blood, respectively .
We have performed quality control steps on raw Illumina reads, followed by alignment using BWA and STAMPY, PCR duplicate marking together with alignment metric calculation by Picard, multi sequence local realignment and base quality score recalibration by GATK as described previously in (2).

Somatic SNP/INDEL and CNV Analysis:
Haplotype caller implemented in Genome Analysis Toolkit (GATK, version 2.5)(3) was used to call variants in tumor--normal pairs. Later we calculated a somatic score according to the method described by Li(4). Variant annotation was performed using Variant Effect Predictor and filtering was performed as detailed in (2).
We performed the CNV analysis on all tumors using the ExomeCNV package. 2

Exome SequenSomatic Structural Variation Analysis:
We have used Breakdancer (5) to call breakpoints. We have used the multi--sample calling mode to assess the somatic status of breakpoints. We have filtered intra-chromosomal breakpoints if: (i) breakpoints were less than 200 bp apart, (ii) there were any supporting reads from the matching blood, (iii) there were less than 8 supporting reads in the tumor, or (iv) breakpoints had a quality score of less than 40. We later annotated the breakpoints using ANNOVAR for hg19 refSeq genes, cyto-band information, segmental duplications, repeat masked regions. Based on these annotations we later filtered breakpoints that overlap with segmental duplications or any repeat elements. We have excluded the breakpoints with supporting reads less than 20 in the Circos plots.

Mutation Signature Analysis:
We have calculated the mutation signature of individual tumors' somatic mutations by considering 6 major mutation classes, i.e., G:C > T:A, G:C > A:T, G:C > C:G, A:T > G:C, A:T > C:G, A:T > T:A. We also calculated the mutation signature using the 5' and 3' based flanking the variation, which leads to 96 classes of signatures.

Clonality Analysis:
Clonality rate is defined to be the percent of tumor cells harboring the identified somatic mutation and correlates with the temporal evolution of the tumor (6, 7).
We estimated the percent of cells that harbor the heterozygous somatic mutations based on the observed variant allele frequency, ploidy at the site of variant and the admixture rate similarly as previously described (8).
We later used the Mclust package in R to cluster the unique somatic mutations scripts. An initial fast alignment was performed on (GRCh37) reference genome followed by local de novo assembly for regions candidate for variation as detailed in (9). Somatic small alterations and CNVs were reported by Complete Genomics (10).
Non--diploid model for the Somatic CNV calls were used from the report, as gliomas can represent gross copy number alterations across the genome. The segments of varying coverage levels were used in Circos plots instead of 100kb merged segments, as focal alterations are critical.
Downstream filtering, integration and analysis scripts were developed in--house.
Somatic structural breakpoints reported in the high confidence dataset by Complete Genomics were further filtered out if (i) frequency in the baseline control set was >