Integrated time course omics analysis distinguishes immediate therapeutic response from acquired resistance

Integrated time course omics analysis distinguishes immediate therapeutic response from acquired... Background: Targeted therapies specifically act by blocking the activity of proteins that are encoded by genes critical for tumorigenesis. However, most cancers acquire resistance and long-term disease remission is rarely observed. Understanding the time course of molecular changes responsible for the development of acquired resistance could enable optimization of patients’ treatment options. Clinically, acquired therapeutic resistance can only be studied at a single time point in resistant tumors. Methods: To determine the dynamics of these molecular changes, we obtained high throughput omics data (RNA-sequencing and DNA methylation) weekly during the development of cetuximab resistance in a head and neck cancer in vitro model. The CoGAPS unsupervised algorithm was used to determine the dynamics of the molecular changes associated with resistance during the time course of resistance development. Results: CoGAPS was used to quantify the evolving transcriptional and epigenetic changes. Applying a PatternMarker statistic to the results from CoGAPS enabled novel heatmap-based visualization of the dynamics in these time course omics data. We demonstrate that transcriptional changes result from immediate therapeutic response or resistance, whereas epigenetic alterations only occur with resistance. Integrated analysis demonstrates delayed onset of changes in DNA methylation relative to transcription, suggesting that resistance is stabilized epigenetically. Conclusions: Genes with epigenetic alterations associated with resistance that have concordant expression changes are hypothesized to stabilize the resistant phenotype. These genes include FGFR1, which was associated with EGFR inhibitors resistance previously. Thus, integrated omics analysis distinguishes the timing of molecular drivers of resistance. This understanding of the time course progression of molecular changes in acquired resistance is important for the development of alternative treatment strategies that would introduce appropriate selection of new drugs to treat cancer before the resistant phenotype develops. Keywords: Acquired resistance, Precision medicine, Data integration, Time course analysis, Genomics, Epigenetics * Correspondence: Christine.Chung@moffitt.org; ejfertig@jhmi.edu Genevieve Stein-O’Brien, Luciane T. Kagohara and Sijia Li contributed equally to this work. Department of Oncology, Sidney Kimmel Comprehensive Cancer Center, Johns Hopkins University, Baltimore, MD, USA Full list of author information is available at the end of the article © The Author(s). 2018 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated. Stein-O’Brien et al. Genome Medicine (2018) 10:37 Page 2 of 22 Background untreated controls, enabling the distinction of the molecu- Recent advances to identification of gene regulation in lar mechanisms associated with acquired resistance from cancer have enabled the selection of targeted therapies those that would occur due to the long-term culturing to inhibit specific regulators of oncogenic signaling path- over the time period that resistance develops. ways essential for tumor development and maintenance Along with novel time course datasets, inferring the [1]. These therapies prolong survival but are not cura- specific and targetable signaling changes that drive tive, since most patients develop acquired resistance therapeutic resistance also requires new bioinformatics within the first few years of treatment [2]. Although a pipelines to analyze and visualize these data. The bio- wide variety of molecular alterations that confer resist- informatics pipelines must integrate genetic, epigenetic, ance to the treatment have been described, the mecha- and transcriptional changes from multiple-high through- nisms and timing of their evolution are still poorly put platforms to infer the complex gene regulatory characterized [3, 4]. As serial biopsies along the treat- mechanisms that are responsible for acquired resistance. ment period are impractical due to the invasiveness and Current supervised bioinformatics algorithms that find high costs of the procedure, the molecular alterations as- time course patterns in genomic data adjust linear sociated with acquired resistance are only known when models to correlate molecular profiles with known tem- resistance has already developed and little is known poral patterns [12–15]. Many unknown variables such as about what changes occur at earlier or later time points culture conditions, immediate response to cetuximab, during the targeted therapy. The lack of adequate in and adaptive changes may have confounding effects on vitro and in vivo time course datasets makes it challen- known covariates of therapeutic response such as growth ging to delineate the two predominant hypotheses for rates, colony size, or apoptosis rates. Unsupervised bio- how therapeutic resistance develops: (1) the presence of informatics algorithms learn the dynamics directly from small populations of resistant cells that will survive the the high-throughput data, and therefore do not require a treatment and repopulate the tumor; or (2) the develop- priori knowledge of the complex dynamics associated ment of de novo resistance mechanisms by the tumor with therapeutic response. Some unsupervised algo- cells [4, 5]. Characterization of the dynamics of genomic rithms [16–21] seek breaking points of coherent, regula- alterations induced during acquired resistance can iden- tory relationships to infer the dynamics along pathways. tify targetable oncogenic drivers and determine the best Many of these algorithms trace individual phenotypes or time point to introduce alternative therapeutic strategies individual genomics platforms. Their ability to determine to avoid resistance establishment [6]. drivers of gene expression associated with acquired re- Epidermal growth factor receptor (EGFR) inhibitors sistance from time course data in multiple experimental represent a common class of targeted therapeutics. conditions and multiple genomics data modalities is Cetuximab, a monoclonal antibody against EGFR, is emerging [22]. Further extensions are needed to contrast FDA approved for the treatment of metastatic colorectal the dynamics of signaling response to therapy to the dy- cancer and head and neck squamous cell carcinoma namics of control conditions to distinguish the specific (HNSCC) [7]. As with other targeted therapies, stable molecular processes that are unique to resistance. response is not observed for a long period and virtually Matrix factorization algorithms are unsupervised and all patients invariably develop acquired resistance [8]. can distinguish the relative molecular changes in each Recent advances in the establishment of in vitro models experimental condition over time without requiring of acquired cetuximab resistance [9] provide a unique prior knowledge of gene regulation. We have found that opportunity to study the time course of genetic events Bayesian, non-negative matrix factorization algorithms resulting in acquired resistance. Cell lines chronically ex- such as Coordinate Gene Activity in Pattern Sets posed to the targeted agent develop resistance and can (CoGAPS) [23] can extend beyond clustering to robustly be sequentially collected during the course of treatment quantify the dynamics and infer the gene regulatory net- to evaluate the progressive molecular changes. Previous works directly from the input time course data [24]. The studies to assess the mechanisms of acquired cetuximab CoGAPS error model can also be modified to enable data- resistance have been limited to comparing the genomic driven inference in distinct molecular platforms for infer- profile of the parental intrinsic sensitive cell line to ence of epigenetic regulation of gene expression [25]. stable clones with acquired resistance [9–11]. Therefore, In this study, we developed a new bioinformatics ana- these studies fail to capture the dynamics of acquired mo- lysis pipeline for integrated analysis of gene expression lecular alterations during the evolution of therapeutic re- and DNA methylation changes that occur during the sistance. The development of in vitro time course data to time course progression of resistance to targeted therap- determine the molecular drivers of therapeutic resistance ies using CoGAPS. Genes uniquely associated with these is crucial. These experimental systems have the further ad- changes were selected using a PatternMarker statistic vantage that time course data can also be generated for [26] to enable novel visualization of molecular Stein-O’Brien et al. Genome Medicine (2018) 10:37 Page 3 of 22 alterations dynamics inferred with CoGAPS. In order to three days for 11 weeks (generations G1 to G11). On the benchmark our new bioinformatics pipeline, we used an eighth day, cells were harvested. Sixty thousand cells in vitro HNSCC cell line model to induce resistance and were replated for another week of treatment with cetuxi- measure the molecular changes using high-throughput mab and the remaining cells were separately collected assays while the resistant phenotype developed. Gene ex- for: (1) RNA isolation (gene expression analysis); (2) pression and DNA methylation changes were screened DNA isolation (DNA methylation analysis); (3) prolifera- weekly while acquired cetuximab resistance was induced tion assay; and (4) storage for future use. All steps were in SCC25 cell line (intrinsic sensitive to cetuximab) and repeated for a total of 11 weeks. In parallel with the compared to the status of the untreated controls at the cetuximab treated cells, we generated controls that re- same culturing time point. CoGAPS [26] inferred spe- ceived the same correspondent volume of phosphate cific patterns of expression and DNA methylation that buffered saline (PBS). Cells were plated in several repli- are associated with the gradual establishment of ac- cates each time at the same initial density. The replicates quired cetuximab resistance. The onset of methylation were then harvested and pooled to provide enough cells changes associated with resistance is temporally delayed for genetic, epigenetic, and proliferation assays. To relative to expression changes suggesting that epigenetic achieve adequate final cell confluence and number of alterations stabilize the transcriptional changes relevant cells for the experimental analysis of each generation, to the resistant phenotype. This analysis found anti- cetuximab- and PBS-treated cells were plated in differ- correlated changes between DNA methylation and gene ent flask sizes. Cells treated with cetuximab were plated expression in FGFR1 during acquired therapeutic resist- in multiple T75 (75cm ) flasks (60,000 cells/flask) that ance. Upregulation of FGFR1 has previously been associ- were combined on the eighth day. PBS-treated cells were ated as a mechanism of acquired cetuximab resistance in plated in a single T175 (175cm ) flask (60,000 cells). HNSCC [27–29]. The identification of a canonical marker This design was selected considering the growth of resistance to EGFR inhibitors in this present study cor- inhibition of the earliest cetuximab generations and to roborates the efficacy of our experimental model and control confluence of the PBS controls at the collection analytical algorithm to identify mechanisms of resistance. time (Additional file 1: Figure S1). To our knowledge, this is the first demonstration of the anti-correlation between FGFR1 methylation and expres- Cell proliferation and colony formation assays sion suggestive of its epigenetic regulation in acquired re- Cell proliferation events were measured using the Click-iT sistance to cetuximab. Thus, this pipeline can identify Plus EdU Flow Cytometry Assay Kit Alexa Fluor 488 Pico- mechanisms of gene regulation in acquired resistance lyl Azide (Life Technologies, Carlsbad, CA, USA) accord- from high-throughput, multi-platform time course data. ing to manufacturer’s instructions. The cetuximab The resulting bioinformatics pipeline is poised to infer the generations were considered resistant when the frequency dynamics of acquired resistance from emerging time of proliferating cells was higher than in the PBS control course data with other cancer types and therapeutics. generations. Proliferation curves were generated using lo- cally weighted polynomial regression (lowess) in R. Methods Anchorage-independent growth assay was used to fur- Cell lines and materials ther confirm the development of resistance. The parental SCC25 cells were purchased from American Type SCC25 and the late G10 resistant cells were treated with Culture Collection (ATCC). Cells were cultured in different concentrations of cetuximab 10 nM, 100 nM, Dulbecco’s Modified Eagle’s medium and Ham’s F12 and 1000 nM. Number of colonies was compared to the medium supplemented with 400 ng/mL hydrocortisone same cells treated with PBS. Colony formation assay in and 10% fetal bovine serum and incubated at 37 °C and Matrigel (BD Biosciences, Franklin Lakes, NJ, USA) was 5% carbon dioxide. The parental cell line SCC25 and the performed as described previously [30]. late cetuximab and PBS generation 10 were authenti- cated using short tandem repeat (STR) analysis kit Stable SCC25 cetuximab resistant single clones (CTXR clones) PowerPlex16HS (Promega, Madison, WI, USA) through Resistance to cetuximab was induced in an independent the Johns Hopkins University Genetic Resources Core passage of SCC25 cells. After resistance was confirmed, Facility. Cetuximab (Lilly, Indianapolis, IN, USA) was single cells were isolated and grown separately to generate purchased from the Johns Hopkins Pharmacy. the isogenic resistant single cell clones (CTXR). In total, 11 CTXR clones were maintained in culture without Induction of cetuximab resistance and time course sample addition of cetuximab. With the exception of one clone collection (CTXR6), all CTXR clones presented substantial survival The HNSCC cell line SCC25 (intrinsically sensitive to advantage compared to the parental SCC25, as reported cetuximab) was treated with 100 nM cetuximab every by Cheng et al. [31]. Each of these clones was Stein-O’Brien et al. Genome Medicine (2018) 10:37 Page 4 of 22 authenticated using STR analysis kit GenePrint 10 HumanMethylation450 BeadChip platform (Illumina) at (Promega) through the JHU-GRCF, as previously the JHMI Sidney Kimmel Cancer Center Microarray published [31]. Core Facility. Briefly, DNA quality was assessed using Proliferation assay was performed to confirm cetuximab the PicoGreen DNA Kit (Life Technologies) and 400 ng resistance in the CTXR clones compared to the parental of genomic DNA was bisulfite converted using the EZ SCC25. A total of 1000 cells were seeded in 96-well plates DNA Methylation Kit (Zymo Research, Irvine, CA, in quadruplicate for each condition. PBS or cetuximab USA) following manufacturer’s recommendations. A (10 nM, 100 nM or 1000 nM) was added after 24 and total volume of 4 μL of bisulfite-converted DNA was 72 h and cells were maintained in culture for seven days. denatured, neutralized, amplified, and fragmented ac- AlamarBlue reagent (Invitrogen, Carlsbad, CA, USA) at a cording to the manufacturer’s instructions. Finally, 12 μL 10% final concentration was incubated for 2 h and fluores- of each sample were hybridized to the array chip cence was measured according to the manufacturer’s rec- followed by primer-extension and staining steps. Chips ommendations (545 nm excitation, 590 nm emission). were image-processed in the Illumina iScan system. Data Resistance in the CTXR clones was confirmed when the from the resulting iDat files were normalized with fun- proliferation rates were higher than in the PBS-treated norm implemented in the R/Bioconductor package minfi SCC25 cells. (version 1.16.1) [35]. Methylation status of each CpG site was computed from the signal intensity in the methyl- RNA-sequencing (RNA-seq) and data normalization ated probe (M) and unmethylated probe (U)asa β value RNA isolation and sequencing were performed for the as follows: parental SCC25 cells (G0) and each of the cetuximab and PBS generations (G1 to G11) and the CTXR clones M β ¼ : at the Johns Hopkins Medical Institutions (JHMI) Deep M þ U Sequencing & Microarray Core Facility. RNA-seq was also performed for two additional technical replicates of Annotations of the 450K probes to the human genome parental SCC25 cell line to distinguish technical variabil- (hg19) were obtained from the R/Bioconductor package ity in the cell line from acquired resistance mechanisms. FDb.InfiniumMethylation.hg19 (version 2.2.0). Probes on Total RNA was isolated from 1 × 10 cells using the sex chromosomes or annotated to single nucleotide AllPrep DNA/RNA Mini Kit (Qiagen, Hilden, Germany) polymorphisms were filtered from analysis. The CpG is- following manufacturer’s instructions. The RNA land probe located closest to the transcription start site concentration was determined by the spectrophotometer was selected for each gene. Genes with CpG island Nanodrop (Thermo Fisher Scientific, Waltham, MA, probes < 200 bp from the transcription start site were USA) and quality was assessed using the 2100 retained to limit analysis to CpG island promoter probes Bioanalyzer (Agilent, Santa Clara, CA, USA) system. An for each gene. Probes were said to be unmethylated for RNA Integrity Number (RIN) of 7.0 was considered as β < 0.1 and methylated for β > 0.3 based upon thresholds the minimum to be used in the subsequent steps for defined in TCGA analyses [34]. All DNA methylation RNA-seq. Library preparation was performed using the data from this study are available from GEO (GSE98813) TrueSeq Stranded Total RNAseq Poly A1 Gold Kit (Illu- as part of SuperSeries GSE98815. mina, San Diego, CA, USA), according to manufacturer’s recommendations, followed by messenger RNA (mRNA) Hierarchical clustering and CoGAPS analysis enrichment using poly(A) enrichment for ribosomal The following filtering criterion for genes from the profil- RNA (rRNA) removal. Sequencing was performed using ing of the time course data from generations of cetuximab the HiSeq platform (Illumina) for 2 × 100 bp sequen- treated cells was used. Genes from RNA-seq data were se- cing. Reads were aligned to hg19 with MapSplice [32] lected if they had log fold change > 1 between any two and gene expression counts were quantified with RSEM time points of the same condition and < 2 between the [33]. Gene counts were upper-quartile normalized and replicate control samples at time zero (5940 genes). CpG log transformed for analysis following the RSEM v2 island promoter probes for each gene were retained if the pipeline used to normalize TCGA RNA-seq data [34]. gene switched from unmethylated (β <0.1) to methylated All RNA-seq data for the cell line mode in this study are (β > 0.3) in any two samples of the time course (1087 available from GEO (GSE98812) as part of SuperSeries genes). We used the union of the sets of genes retained GSE98815. from these filtering criteria on either data platform for analysis, leaving a total of 6445 genes in RNA-seq and DNA methylation hybridization array and normalization 4703 in DNA methylation. Genome-wide DNA methylation analysis was performed Hierarchical clustering analysis was performed with on the same samples as RNA-seq using the Infinium Pearson correlation dissimilarities between genes and Stein-O’Brien et al. Genome Medicine (2018) 10:37 Page 5 of 22 samples on all retained genes. CoGAPS analysis was per- meta-pathway were derived for each pattern using the formed on both log transformed RNA-seq data and PatternMarkers statistics [26]. Comparisons between DNA methylation β values, independently using the R/ DNA methylation and gene expression values for Bioconductor package CoGAPS [23] (version 2.9.2). PatternMarkerGenes or from CoGAPS patterns and am- CoGAPS decomposes a matrix of data D according to plitudes were computed with Pearson correlation. the model Gene set analysis of cetuximab resistance signatures, the EGFR network, and pathways D ∼N A P ; Σ ; i; j i;k k; j i; j Gene set activity was estimated with the gene set statis- k¼1 tic implemented in calcCoGAPSStat of the CoGAPS R/ Bioconductor package [23]. Analyses were performed on where N represents a univariate normal distribution, three gene sets: resistance signatures, gene targets of matrices A and P are learned from the data for a speci- transcription factors in the EGFR network, and canon- fied number of dimensions P, Σ is an estimate of the i, j ical pathways. Resistance signatures were defined based standard deviation of each row and column of the data on previous literature. Specifically, in a previous study, matrix D, and i represents each gene and j each sample. CoGAPS learned a meta-pathway from gene expression In this decomposition, each row of the pattern matrix P Val12D data corresponding to overexpression of the HRAS quantifies the relative association of each sample with a in the HaCaT cell line model. That study associated the continuous vector of relative gene expression changes in CoGAPS HaCaT-HRAS meta-pathway with gene expres- the corresponding column of A. That is each row of P sion changes in acquired cetuximab resistance in the provides the relative magnitude across samples are called HNSCC cell line UMSCC1 [23]. In the current study, we patterns and quantify the separation of distinct experi- applied the PatternMarkers statistic [26] to the previ- mental conditions. These relative gene weights in the ously published CoGAPS analysis of these data to derive columns of A represent the degree to which each gene is a gene set from this meta-pathway called HACAT_ associated with an inferred pattern and are called meta- HRAS_CETUXIMAB_RESISTANCE or HACAT_RESIS pathways. Together, these matrices provide a low- TANCE. In addition, we searched MSigDB [37] (version dimensional representation that reconstructs the signal 5.2) for all gene sets associated with resistance to EGFR of the input genomics data. A single gene may have inhibition. In this search, we found the gene sets COL non-zero magnitude in several distinct gene sets, repre- DREN_GEFITINIB_RESISTANCE_DN and COLDREN_ senting the fact that a single gene can have distinct roles GEFITINIB_RESISTANCE_UP representing resistance in different biological processes (such as immediate to the EGFR inhibitor gefitinib in non-small-cell lung can- therapeutic response and acquired resistance). A recently cer (NSCLC) cell lines [38]. Gene sets of transcription fac- developed PatternMarker statistic [26] selects the genes tor targets were obtained from experimentally validated that are unique to each of the inferred patterns and targets annotated in the TRANSFAC [39] professional therefore represent biomarkers unique to the corre- database (version 2014.1). Canonical pathways were ob- sponding biological process. tained from the C2 set of MSigDB [37](version6.1). In the CoGAPS analysis of the data in this study, the standard deviation of the expression data was 10% of the Sources and analysis of additional in vitro and human signal with a minimum of 0.5. The standard deviation of tumor genomics data DNA methylation data under the assumption that β Genomics analyses of TCGA were performed on level 3 values follow a beta distribution is RNA-seq and DNA methylation data from the 243 HPV- vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi negative HNSCC samples from the freeze set for publica- β 1−β i; j i; j tion [36]. DNA methylation data were analyzed for the Σ ¼ : i; j M þ U þ 1 same CpG island promoter probes obtained in the cell line i; j i; j studies. Pearson correlation coefficients were computed in CoGAPS was run for a range of 2–10 dimensions P R to associate different molecular profiles. for expression and 2–5 for DNA methylation. Robust- Additional analysis was performed on Affymetrix Hu- ness analysis with ClutrFree [36] determined that the man Genome U133 plus 2.0 GeneChip arrays for the optimal number of dimensions P for expression was 5. SCC1/1CC8 isogenic cetuximab sensitive and resistant DNA methylation was run in four parallel sets using cell line pair described previously (GEO GSE21483 [30]). GWCoGAPS [26]. In DNA methylation, the maximum Additional gene expression data from SCC25 generated number of patterns that modeled resistance mechanisms from the same platform in the same lab was also used over and above technical variation in replicate samples for analysis, using fRMA for normalization [40] to con- of SCC25 was three. Gene sets representative of the trol for batch effects as described previously [41]. Stein-O’Brien et al. Genome Medicine (2018) 10:37 Page 6 of 22 Analysis was also performed on gene expression data 2000). It was prospectively planned to perform transla- measured with Illumina HumanHT-12 WG-DASL V4.0 tional research and patients gave their informed consent R2 expression beadchip arrays on pre-treatment samples for repeated biopsies. from patients subsequently treated with cetuximab from Bossi et al. [42], using expression normalization and Results progression-free survival groups as described in the Prolonged exposure to cetuximab induces resistance study. Data were obtained from the GEO GSE65021 Cetuximab resistance was induced by treating the series matrix file. SCC25 cells for a period of 11 weeks (CTX-G1 to – DNA samples from eight human tumor surgical speci- G11). SCC25 cells treated with PBS were used as time-- men post cetuximab treatment from the sample cohort matched controls (PBS-G1 to –G11). Response to cetux- in Schmitz et al. [43] were obtained for methylation pro- imab was determined by comparing the proliferation filing. Specifically, for each tumor one FFPE slide was rates between CTX and PBS generations. Proliferation of stained with hematoxylin and eosin and tumor burden the PBS generations is stable throughout the 11 weeks was evaluated. When the tumor content was < 50%, the (G1 to G11, Fig. 1a). Conversely, proliferation of the adjacent unstained FFPE slides were macrodissected in CTX generations progressively increases over each week order to enrich the tumor burden. A double code was (Fig. 1a). Relative to the untreated controls, the growth assigned to each sample by Biorepository. DNA was then of the treated cells is initially (CTX-G1) inhibited until extracted from two unstained slides using the QIAamp CTX-G3. Starting at CTX-G4, the absolute proliferation DNA FFPE Tissue kit. Briefly, slides were dipped into a values are equal at this week, but the fit to the data xylene bath until paraffin was melted. Then slides were across all time points suggests that the cells become re- washed with ethanol 100% and tissue was harvested for sistant to the anti-proliferative effects of cetuximab and extraction with QIAGEN affinity columns. The extracted gain stable growth advantages compared to the un- DNA was quantified by NanoDrop spectrophotometer. treated controls (CTX-G8 to –G11). DNA methylation was measured with the Illumina Comparison of proliferation rates between generations MethylationEPIC BeadChip (850K) array. Array data of CTX-treated cells relative to generations of cells were normalized with the NOOB method [44] and con- treated with PBS enabled us to conclude that cell growth verted to virtual 450K arrays using the R/Bioconductor advantages arise from chronic cetuximab treatment and package minifi version 1.22.1 and are available from are associated with resistance rather than prolonged cell GEO SuperSeries (GSE110995). Two samples had DNA culturing. We mirrored the changes in proliferation rates content < 250 ng and clustered separately from the with clinical responses seen in HNSCC tumors treated remaining samples. These samples were filtered as low with cetuximab (Fig. 1b). The lower growth rates in quality and excluded from the analysis, leaving six total CTX-G1 to -G3 may be an equivalent to the initial ef- tumor samples with DNA methylation data. Probes se- fects of the clinical treatment when the cancer cells are lected for the in vitro Illumina 450K DNA methylation sensitive to cetuximab and reduction of the tumor size is data were used for subsequent analyses. Gene expression observable. Even with gain in cell proliferation at CTX- data from biopsy samples before cetuximab treatment G3 and -G4, our model still corresponds to response to and surgical samples after cetuximab treatment were ob- the treatment since the treated generations are not tained from the previous Schmitz et al. [43] study and growing more than the controls (clinical stable tumor normalized as described previously [41] and available size). Finally, from CTX-G4 the higher proliferation even from GEO SuperSeries (GSE110996). with cetuximab treatment is a representation of acquired We performed t-tests and projections in R on the resistance noticeable in the HNSCC patients as tumor probe that had the highest standard deviation of expres- recurrence or increase in tumor size. sion values for each gene. CoGAPS signatures were also Higher proliferation in treated than in untreated cells projected into these gene expression data using the starts at CTX-G4 and we established this time point to methods described in Fertig et al. [23] with the ProjectR call as the moment at which cetuximab resistance is sta- package version 0.99.15 available from Github (https:// bly acquired and all subsequent time points continue to github.com/genesofeve/projectR). We also performed t- develop acquired stable cetuximab resistance. To con- tests in R to compare the long- (LPFS) and short-term firm this hypothesis, we evaluated the ability of the re- progression-free survival (SPFS) groups based on the sistant CTX-G10 to anchorage-independent growth. values obtained from this projection. Even under different concentrations of cetuximab, HNSCC samples and patient information collection CTX-G10 presents enhanced anchorage-independent were approved by the Independent Ethics Committee growth compared to the parental SCC25 (G0) (two-way and the Belgian Health Authorities and conducted in ac- ANOVA with multiple comparisons p value < 0.01 for cordance with the Declaration of Helsinki (October each concentration, Additional file 1: Figure S2), Stein-O’Brien et al. Genome Medicine (2018) 10:37 Page 7 of 22 a Proliferation assay Generations (weeks of treatment) b Clinical evolution Tumor detection Treatment with cetuximab Tumor reduction Stable tumor size Tumor growth / Disease recurrence RESPONSE TO CETUXIMAB ACQUIRED RESISTANCE Unaccessible molecular changes Fig. 1 In vitro time course reflects clinical evolution of cetuximab response and evolution of acquired resistance. Intrinsic cetuximab-sensitive HNSCC cell line SCC25 was treated with cetuximab (red) or PBS (black) for 11 generations to develop acquired resistance. Proliferation assay (flow) of cetuximab treatment (red line) and PBS-treated cells (black line) measured cetuximab response for all SCC25 generations. While proliferation of the PBS generations was stable throughout the 11 weeks, proliferation of the CTX generations progressively increased over each week. Relative to the untreated controls, the growth of the treated cells was initially (CTX-G1) inhibited until CTX-G3. Starting at CTX-G4, the cells became resistant to the anti-proliferative effects of cetuximab and gained stable growth advantages compared to the untreated controls demonstrating the stabilization of cetuximab resistance proliferation rates than the PBS controls (shown in in later generations. Fig. 1a). The two groups of samples resistant to cetuximab are represented by a progressive increase Treatment vs control gene expression changes governs in proliferation that is more significant than in the clustering and immediate therapeutic response is untreated controls (weeks 4 to 8) and by the stabilization confounded with changes from acquired resistance in the proliferation rates (weeks 9 to 11), but still higher RNA-seq data for the parental SCC25 cell line (G0) and than in thePBS generations. Theexpressionchanges at from each generation of CTX- and PBS-treated cells the distinct time points during development of acquired were collected to characterize the gene expression resistance are shared among numerous genes. Although changes occurring as cells acquired cetuximab resist- the clustering was able to separate cetuximab- from ance. Gene expression changes between treated (cetuxi- PBS-treated cells, it was not possible to discriminate mab) and untreated (PBS) cells and over generations of the alterations related to an immediate therapeutic re- treated cells are apparent in time-ordered RNA-seq data sponse (not relevant to the resistant phenotype) from (Fig. 2a). Additional clustering analysis of the samples resistance-specific gene expression changes. accounting for the treatment time point (generations/ Similar separation of the sensitive and resistant phases columns) (Additional file 1: Figure S3) distinguish three of cetuximab response is observed in clustering analysis clusters of samples: those with cetuximab sensitivity of gene signatures previously described in HNSCC and (CTX-G1 to CTX-G3); those with early cetuximab re- NSCLC cell line models resistant to cetuximab or gefi- sistance (CTX-G4 to CTX-G8); and those with late or tinib (anti-EGFR small molecule), respectively [38, 41] stable cetuximab resistance (CTX-G9 to CTX-G11). The (Additional file 1: Figure S4). For these genes, changes group of cetuximab sensitive samples corresponds to the during early phases of resistance clusters for CTX-G4 to time points at which the CTX generations present lower CTX-G6 as distinct from later generations CTX-G7 to Cell proliferation (%) 10 15 20 25 30 35 40 Stein-O’Brien et al. Genome Medicine (2018) 10:37 Page 8 of 22 Clustering of gene expression data CoGAPS analysis of gene expression data PBS Cetuximab Expression of PatternMarker genes bc CoGAPS expression patterns for CoGAPS expression patterns (sample weights v. time) PBS Cetuximab 0 24 68 10 Color Key Row Z-score (log2 RSEM) Weeks of treatment PBS −4 0 4 Cetuximab weeks of treatment Fig. 2 CoGAPS analysis identifies signatures of resistance to EGFR inhibitors and separate resistant and control generations. a Heatmap of gene expression values in 11 generations of SCC25 cells treated with 100 nM of cetuximab (red columns) to acquire resistance and with PBS as control (black columns). b Heatmap of gene expression values for PatternMarker genes identified with CoGAPS analysis of gene expression data from 11 generations of SCC25 cells treated with PBS as control (black columns) and with 100 nM of cetuximab (red columns) to acquire resistance. Rows are colored according to which CoGAPS pattern the PatternMarker statistic assigned each gene and sorted by the PatternMarker statistic. c CoGAPS patterns inferred from gene expression data over generations of PBS control (black lines) or treatment with 100 nM of cetuximab (red lines) CTX-G11. Nevertheless, these signatures also cluster gene can have high amplitude in multiple patterns, mod- samples with gene expression changes at early phases eling multiple regulation of genes but complicating (CTX-G1 to CTX-G3) as distinct from samples from visualization of the inferred patterns. A recently devel- PBS-treated generations. However, these analyses were oped PatternMarker statistic defines genes that are insufficient to quantify the relative dynamics of genes as- uniquely associated with each of these patterns. Limiting sociated with immediate response to therapy or subse- the heatmap to these genes enables visualization of the quent acquired resistance. dynamics of gene expression changes in our time course dataset (Additional file 1: Figure S5). CoGAPS analysis of gene expression distinguishes patterns In this heatmap of CoGAPS PatternMarker genes, we of acquired resistance from immediate therapeutic response observe that only three patterns (EP1, EP2, and EP3) dis- To define gene expression signatures for treatment effect tinguish the experimental conditions (cetuximab vs PBS) and cetuximab resistance, we applied the CoGAPS [26] (top three patterns on Additional file 1: Figure S5). The Bayesian matrix factorization algorithm to the time other two patterns, EP4 and EP5, represent changes in course gene expression data. CoGAPS decomposes the gene expression from the parental cell lines and subse- input data into two matrices: a pattern matrix with rela- quent generations or an expression pattern that is con- tive sample weights along rows and an amplitude matrix stant and corresponds to signature of highly expressed with relative gene weights along columns. Each row of genes (lower two patterns on Additional file 1: Figure the pattern matrix quantifies the extent of transcrip- S5), respectively. We note that highly expressed genes tional changes within the genes in the corresponding associated with EP5 may also have dynamic changes due column of the amplitude matrix and provides a low di- to treatment and are filtered in the PatternMarker ana- mensional representation of the biological process in lysis of all patterns in Additional file 1: Figure S5. EP4 that data. In this analysis, we identified five CoGAPS represents expression changes between treated cells and patterns (Expression Patterns [EP]) (Additional file 1: the parental cell line, which have a technical effect on Figure S5) in the time course gene expression dataset. A gene expression. Notably, even the exclusion the flat 11 Stein-O’Brien et al. Genome Medicine (2018) 10:37 Page 9 of 22 pattern for highly expressed genes (EP5) still retains expression signatures associated consistently with the highly expressed genes in the signature. Retaining the development of cetuximab resistance. Gene set statistics technical pattern (EP4) in this calculation filters genes of transcription factor targets of EGFR on CoGAPS gene with expression changes from technical artifacts in the weights are significantly downregulated in this acquired experimental conditions. The resulting set of Pattern- resistance pattern (Additional file 1: Figure S6). One strik- Marker genes for EP1–EP3 enable visualization of the ing exception is c-Myc, which trends with acquired resist- expression changes dynamics that are associated with ance (p value of 0.06), consistent with the role of this cetuximab response (Fig. 2b) and allow the definition of transcription factor in cellular growth. Resistance signa- a gene signature associated with that response ture COLDREN_GEFITINIB_RESISTANCE_DN is sig- (Additional file 2: Table S1). nificantly downregulated in EP2 (p value of 0.04). There Similar to the separation seen with clustering are 32 statistically significant canonical pathways associ- (Additional file 1: Figure S5), the CoGAPS pattern EP1 ated with this pattern, including notably telomerase, PI3K, distinguishes cetuximab from PBS at every generation and cell cycle pathways (Additional file 3: Table S2). (Fig. 2b and c, top). These genes present an immediate The third CoGAPS expression pattern (EP3) repre- transcriptional upregulation in response to cetuximab sents a gradual repression of gene expression with treatment. Gene set analysis to determine the function cetuximab treatment (Fig. 2b and c, bottom). This of CoGAPS patterns was performed with an enrichment expression pattern trends to significant enrichment in the analysis on all gene weights in the amplitude matrix ob- COLDREN_GEFITINIB_RESISTANCE_DN resistance tained from the CoGAPS analysis. By performing the signature (Additional file 1:Figure S6,one-sided p value 0. analysis on gene weights and not only the PatternMarker 12) and downregulated in the HACAT_HRAS_CETUXI- genes, as shown in Fig. 2b, we account for multiple regu- MAB_RESISTANCE resistance signature (Additional file 1: lation of genes in pathways. Specifically, we performed Figure S6, one-sided p value 0.09). This confirms that EP3 gene set analysis on published resistance signatures [38, is associated with repression of gene expression during ac- 41], transcription factors previously associated with the quired cetuximab resistance. There are also 29 statistically EGFR signaling network during cetuximab response in significant canonical pathways associated with this pat- HNSCC [39, 41], and canonical pathways from MSigDB tern, including cell lineage, metabolic, WNT, and GSK3 [37, 45] (Additional file 1: Figure S6; Additional file 3: pathways (Additional file 3: Table S2). Table S2). Gene set analysis confirms that published re- sistance signatures [38, 41] are significantly enriched in EP1 (Additional file 1: Figure S6; one-sided p values of 0. Changes in DNA methylation inferred with CoGAPS are 002 and 0.003 for resistance gene sets COLDREN_GEFI- associated with resistance to cetuximab, but not the TINIB_RESISTANCE_DN and HACAT_HRAS_CETUX immediate response to treatment observed in gene IMAB_RESISTANCE, respectively). However, the tran- expression scriptional changes in this pattern are not associated To determine the timing of the methylation changes as- with acquired resistance to cetuximab and even decrease sociated with acquired resistance, we also measured modestly as resistance developed. Further, enrichment DNA methylation in each cetuximab generation of by transcription factor AP-2alpha targets (TFAP2A; one- SCC25 cells and PBS controls (Fig. 3a). Application of sided p value of 0.05) confirms previous work indicating the CoGAPS matrix factorization algorithm with the that transcription by AP-2alpha is induced as an early PatternMarker statistics to the methylation data reveals feedback response to EGFR inhibition [39]. There are 84 a total of three methylation patterns (MP) (Fig. 3bc; significant canonical pathways from MSigDB, including Additional file 2: Table S1): gradual increase of DNA notably pathways associated with the immune system, methylation in controls (MP1, Fig. 3b middle); rapid extracellular matrix, ERBB4 signaling, and VEGF signal- demethylation in CTX generations starting at CTX-G4 ing (Additional file 3: Table S2). Based upon these (MP2, Fig. 3b bottom); and rapid increase in DNA findings, we concluded that EP1 is associated with im- methylation in CTX generations starting at CTX-G4 mediate response to cetuximab although it includes (MP3, Fig. 3b top). In contrast to the gene expression genes that are also associated with cetuximab resistance data, there is no immediate shift in DNA methylation in previous studies. resulting from cetuximab treatment. Gene set analysis The second CoGAPS expression pattern (EP2) quanti- was performed on canonical pathways from MSigDB fies divergence of the cetuximab treated cells from (Additional file 4: Table S3) and found 26 statistically controls at generation CTX-G4 (Fig. 2b and c, middle) significant pathways for MP1, 29 for MP2, and 27 for which is the time point that cetuximab treated cells MP3. In contrast to gene expression, the majority of ca- present significant and stable growth advantage over nonical pathways is shared by the three methylation PBS controls (Fig. 1a). Therefore, EP2 contains gene patterns and include notably the cytokine (PID-IL8- Stein-O’Brien et al. Genome Medicine (2018) 10:37 Page 10 of 22 Clustering of DNA methylation data CoGAPS analysis of DNA methylation data PBS Cetuximab DNA methylation of PatternMarker genes CoGAPS DNA methylation patterns bc for CoGAPS methylation patterns (sample weights v. time) PBS Cetuximab 1 1 weeks of treatment 0 24 68 10 Weeks of treatment Color Key Row Z-score (b) PBS Cetuximab −4 0 4 weeks of treatment Fig. 3 Dynamics of DNA methylation alterations and patterns in acquired cetuximab resistance. a Heatmap of DNA methylation values in 11 generations of SCC25 cells treated with PBS as control (black columns) and with 100 nM of cetuximab (red columns) to acquire resistance. b Heatmap of DNA methylation values for genes selected by CoGAPS DNA methylation patterns analysis in the same SCC25 cetuximab and PBS generations. c CoGAPS patterns inferred from DNA methylation data over generations of PBS control (black lines) or treatment with 100 nM of cetuximab (red lines) CXCR2 and IL8-CXCR1 pathways) and FGFR of the expression and methylation patterns previously (Reactome Signaling by FGFR3 mutants) signaling shown (Figs. 2c and 3c, respectively) enable visualization pathways. of the time point when changes between cetuximab and Comparing the CoGAPS patterns from gene expression PBS generations are significant in each pattern. These and DNA methylation reveals strong anti-correlation be- dynamics explain the discrepancy between the genes as- tween gene expression and DNA methylation in resistant sociated with each pattern and suggest that DNA methy- patterns (Additional file 1: Figure S7A). We observed that lation stabilizes the gene expression signatures crucial to the gene expression changes associated with acquired re- the maintenance of acquired cetuximab resistance. sistance occur gradually and are evident in early genera- tions (Fig. 2c). The DNA methylation is consistent in Gene expression and methylation profile of SCC25 cetuximab treatment and control PBS in DNA methyla- single-cell clones with acquired cetuximab resistance tion patterns MP2 and MP3 during early generations; fol- demonstrates cell heterogeneity lowing which, there is a rapid accumulation in DNA The little overlap between the gene expression and DNA methylation changes starting after generations CTX-G4 methylation PatternMarker genes and non-specific DNA and CTX-G5 in both MP2 and MP3 (Fig. 3c), concurrent methylation pathways may arise due to the development with the onset of the observed growth advantage over the of different resistant sub-clones with specific gene signa- PBS control (Fig. 1a). tures of acquired resistance in bulk data. In order to While the patterns themselves are anti-correlated, the address this issue and to delineate whether our pre- gene weights that define meta-pathways and are inferred sumptive drivers resulted from clonal expansion of re- in the amplitude matrix corresponding to each pattern sistant cells or from the development of new epigenetic with CoGAPS are not (Additional file 1: Figure S7B). alterations to drive resistance, we measured DNA We also observed little overlap between the PatternMar- methylation and gene expression on a panel of 11 iso- ker genes from methylation patterns and gene expres- genic stable cetuximab-resistant clones (CTXR1 to sion. Changes in DNA methylation are delayed relative CTXR11) derived from SCC25 cells in a previous study to those of gene expression in acquired cetuximab resist- [31]. Despite being derived from the parental SCC25 ance as can be noted in Fig. 4, where direct comparison cells after chronic exposure to cetuximab, the CTXR 11 Stein-O’Brien et al. Genome Medicine (2018) 10:37 Page 11 of 22 Fig. 4 DNA methylation and expression CoGAPS patterns demonstrate delayed onset of epigenetic changes in acquired resistance. CoGAPS patterns for gene expression (top) and DNA methylation (bottom) of patterns associated with acquired cetuximab resistance in SCC25 cetuximab generations (red) relative to PBS generations (black). Vertical dashed line represents time at which patterns for SCC25 generation separated from pattern for PBS generations. The timing of methylation changes distinguishing cetuximab-resistant generations was delayed in DNA methylation relative to that of gene expression clones and the time course generations display wide- analysis between DNA methylation and gene expression for spread differences. We plot both expression and methy- each of the DNA methylation PatternMarker genes (Fig. 5). lation profiles (Additional file 1: Figure S8 and Figure This analysis identified FGFR1 as one of the genes S9, respectively) among the DNA methylation Pattern- with significant anti-correlation between expression Marker genes that are anti-correlated with expression and methylation, suggesting potential epigenetic regu- and cellular morphology (Additional file 1: Figure S10) lation during cetuximab resistance acquisition. This for these clones. Significantly greater heterogeneity is finding is consistent with previous studies that associ- observed among the CTXR clones in all these genomic ate differential expression of FGFR1 with resistance to and morphological data. Figure 5 demonstrates that EGFR inhibitors, including cetuximab, in HNSCC and higher heterogeneity among single cell clones is also ob- other tumor types in vitro and in vivo [27, 46–48]. served in the epigenetically regulated PatternMarker However, none of these studies demonstrate an asso- genes from the CoGAPS analysis. These results suggest ciation between FGFR1 upregulation and demethyla- that different mechanisms of resistance may arise in the tion. Given the tight temporal regulation of the DNA same HNSCC cell line as a result of intra-heterogeneity, methylation PatternMarkers with anti-correlated expres- resulting in the detection of a wide range of expression sion and the previous work on FGFR1,we hypothesize signatures with higher or lower correlation with the that this set of genes represents epigenetic drivers of methylation profile depending on the size of each spe- acquired resistance. cific cell population. We hypothesize that epigenetically regulated genes shared along the time course patterns and resistant FGFR1 overexpression and demethylation are associated single-cell clones might implicate common mechanisms with acquired cetuximab resistance in the time course acquired during evolution of the stable resistance pheno- and in stable cetuximab-resistant clones type. To test this assumption, we also performed correl- To ascertain potential drivers of the stable cetuximab- ation analysis for each of the epigenetically regulated resistant phenotype induced by DNA methylation, we genes in our resistant set (Fig. 5) in the resistant single defined genes that are PatternMarkers [26] of the DNA cell clones and parental cell lines. Nine of the epigeneti- methylation patterns associated with stable acquired cally regulated PatternMarker genes also have significantly cetuximab resistance (MP2 and MP3). We then applied anti-correlated gene expression and DNA methylation in correlation analysis to determine genes that were epi- the stable cetuximab-resistant clones. Of these, only genetically regulated. Specifically, we performed correlation FGFR1 is demethylated and re-expressed in a CTXR Stein-O’Brien et al. Genome Medicine (2018) 10:37 Page 12 of 22 PatternMarker genes for CoGAPS methylation patterns that are anti-correlated with expression in time course data Gene expression in time course and resistant clones DNA methylation in time course and resistant clones ab PBS Cetuximab Resistant clones PBS Cetuximab Resistant clones EVC2 EVC2 PIGG PIGG ZNF721 ZNF721 RRAD RRAD FRYL FRYL ZNF141 ZNF141 GALC GALC GNPDA2 GNPDA2 MAEA MAEA MAEA MAEA ZNF229 ZNF229 LMX1B LMX1B VAMP5 VAMP5 ZNF71 ZNF71 PRIMA1 PRIMA1 IDUA IDUA ZNF569 ZNF569 KCNJ5 KCNJ5 GSTA4 GSTA4 ITGA6 ITGA6 HOPX HOPX HOPX HOPX KLF3 KLF3 SLC39A8 SLC39A8 SLC39A8 SLC39A8 HLA−DPB1 HLA−DPB1 CFLAR CFLAR ENPEP ENPEP SOCS6 SOCS6 SLC35D1 SLC35D1 LCORL LCORL SLC2A10 SLC2A10 NPNT NPNT TLL1 TLL1 RDM1 RDM1 NEUROG2 NEUROG2 ECHDC2 ECHDC2 GAB1 GAB1 ETS2 ETS2 DOCK10 DOCK10 GATM GATM RIPK3 RIPK3 HEYL HEYL SLC6A2 SLC6A2 COL2A1 COL2A1 ASRGL1 ASRGL1 MTTP MTTP MIR155HG MIR155HG TFAP2E TFAP2E ZNF483 ZNF483 ALS2CL ALS2CL GNAL GNAL FBXO15 FBXO15 DIO3 DIO3 BACH1 BACH1 SYT5 SYT5 LRP10 LRP10 ELOVL4 ELOVL4 COL24A1 COL24A1 THSD1 THSD1 KCNE3 KCNE3 FLT4 FLT4 ELMOD3 ELMOD3 TSPAN10 TSPAN10 SDCBP2 SDCBP2 TTC21A TTC21A RYR3 RYR3 ISLR2 ISLR2 CNPY1 CNPY1 WNT4 WNT4 EFCAB10 EFCAB10 IGFBPL1 IGFBPL1 IL17RB IL17RB GDF15 GDF15 GAS7 GAS7 LTBP3 LTBP3 RAD50 RAD50 KCNJ2 KCNJ2 PPARG PPARG SOX15 SOX15 FGFR1 FGFR1 HS6ST1 HS6ST1 KIFC3 KIFC3 CASP3 CASP3 IQGAP2 IQGAP2 ODC1 ODC1 ODC1 ODC1 C2orf70 C2orf70 ICAM5 ICAM5 PFKP PFKP PFKP PFKP PYCARD PYCARD CIT CIT STC2 STC2 TOX TOX TCFL5 TCFL5 FBN2 FBN2 FBN2 FBN2 FOXK2 FOXK2 SPSB4 SPSB4 GAPDHS GAPDHS LBX1 LBX1 TMEM145 TMEM145 COL8A1 COL8A1 GRB10 GRB10 CAD CAD ABCF2 ABCF2 WBP2NL WBP2NL TMEM45B TMEM45B BCR BCR TIMM50 TIMM50 Row Z-score (β) Row Z-score Weeks of treatment Weeks of treatment −4 0 4 −4 0 4 Fig. 5 Clonal heterogeneity does not reflect signature of epigenetically regulated genes observed in bulk time course analysis. a Heatmap of gene expression values for DNA methylation PatternMarker genes for acquired resistance that were anti-correlated between expression and DNA methylation (Fig. 4). Data includes 11 generations of SCC25 cells treated with PBS as control (black columns labeled PBS) and with 100 nM of cetuximab (red columns labeled cetuximab) to acquire resistance and gene expression data from independent, stable cetuximab-resistant clones in absence of cetuximab treatment (CTX-resistant clones). Gene expression heatmap on a red-yellow scale indicated in the color key. b Heatmap of DNA methylation data in conditions described in (a), on a blue-yellow scale indicated in the color key clone relative to the parental SCC25 cell line (Fig. 6a, pre and post cetuximab treatment, we further investigate remaining eight genes in Additional file 1:FigureS11). the pattern of expression and methylation of FGFR1 and In this analysis, overexpression and de-methylation of EGFR in publicly available datasets. Using gene expres- FGFR1 expression occurs in only one of the resistant sion and DNA methylation data from The Cancer clones (CTXR10). QRT-PCR of FGFR1 gene expression Genome Atlas (TCGA) for 243 HPV-negative HNSCC in CTXR10 relative to the parental cell line demon- pre-treatment samples [36], we verified that the upregu- strated a > 30-fold change (Fig. 6b). Furthermore, in the lation of EGFR and FGFR1 is not concomitant (Pearson resistant cell clone, increased levels of FGFR1 were as- correlation coefficient = − 0.06, p value = 0.33, Fig. 7a). sociated with increased levels of phospho-FGFR1 and Additionally, the negative correlation of FGFR1 gene decrease in EGFR and phospho-EGFR as assessed by expression and DNA methylation status is statistically Western blot (Fig. 6c). This clone is one of the fastest significant (Pearson correlation r of − 0.32, p value growing under cetuximab treatment (Additional file 1: < 0.0001, Fig. 7b),suggesting that FGFR1 transcription is Figure S12). This observation suggests that the bulk associated with demethylation in some HPV-negative data from the time course captured clonal outgrowth of HNSCC tumors. Since there is no treatment information a cetuximab-resistant clone with similar molecular fea- available for the TCGA dataset, we could not make as- tures (FGFR1 demethylation) to CTXR10 and that sumptions related to cetuximab resistance and whether clonal outgrowth is the dominant mechanism of resist- FGFR1 methylation is a consequence of the treatment. To ance in our model. assess this question, we collected new DNA methylation data for six HNSCC tumors after cetuximab treatment Observed FGFR1 dynamics in vitro recapitulates relationships from a cohort of HNSCC tumor samples described previ- from in vivo tumor genomicsand acquired cetuximab ously [43]. All six samples have low DNA methylation resistance values for FGFR1 (β-values in the range of 0.04–0.08, with In order to confirm that the mechanisms we found with a mean of 0.05), suggesting that the gene is unmethylated our in vitro approach are present in HNSCC samples in these samples. While there was insufficient DNA to SCC25 CTXR10 CTXR3 CTXR1 CTXR12 CTXR7 CTXR8 CTXR2 CTXR4 CTXR9 CTXR11 CTXR5 SCC25 CTXR1 CTXR12 CTXR4 CTXR11 CTXR9 CTXR5 CTXR7 CTXR8 CTXR3 CTXR10 CTXR2 Stein-O’Brien et al. Genome Medicine (2018) 10:37 Page 13 of 22 a FGFR1 gene expression vs DNA methylation in resistant clones CTXR10 CTXR1 CTXR12 CTXR11 CTXR3 CTXR4 CTXR2 CTXR8 CTXR9 CTXR7 7 CTXR5 SCC25 0.1 0.2 0.3 0.4 FGFR1 DNA methylation (beta cg20148210) b FGFR1 gene expression c Protein expression SCC25 CTXR10 40 FGFR1 pFGFR (Tyr653/654) 20 EGFR pEGFR (Tyr1068) 0 GAPDH SCC25 CTXR10 Fig. 6 Overexpression and de-methylation of FGFR1 in acquired cetuximab resistance is confirmed in stable SCC25-resistant clones. a Expression of FGFR1 relative to DNA methylation in stable cetuximab-resistant clones. b QRT-PCR of FGFR1 gene expression in CTXR10 relative to the parental cell line (> 30-fold change). c Western blot comparing FGFR1, phospho-FGFR1, EGFR, and phospho-EGFR in CTXR10 relative to the parental SCC25 cell line. In the resistant cell clone, increased levels of FGFR1 were associated with increased levels of phospho-FGFR1 and decrease in EGFR and phospho-EGFR quantify DNA methylation before treatment in these pa- Nonetheless, these findings suggest that similar molecular tients, FGFR1 gene expression increases after treatment in mechanisms may contribute to both mechanisms (intrin- four of the six tumor samples (Additional file 1:Figure sic and acquired) of cetuximab resistance in HNSCC. S13). While this cohort is small, these data and TCGA suggest that FGFR1 methylation is potentially associated CoGAPS signatures of resistance and therapeutic response with its re-expression in HNSCC tumor samples in re- replicate in an independent in vitro system and sponse to cetuximab treatment. significantly stratified patient samples with long- To determine whether FGFR1 is associated with cetux- vs short-progression-free survival imab response, we used gene expression data from To further illustrate that the results are reflective of HNSCC patients before cetuximab treatment available HNSCC in a general fashion, we evaluated the behavior from Bossi et al. [42]. After follow-up, the patients were of the two additional cell lines and human tumors in the separated in SPFS (median three months survival) or CoGAPS signatures using gene expression data available LPFS (median 19 months survival) according to time of from previously published studies. The HNSCC cell lines recurrence or metastasis development. Using this data- SCC1 and 1CC8 were chosen as the cetuximab-resistant set, we confirmed that EGFR expression in SPFS is 1CC8 was generated from the cetuximab sensitive SCC1 significantly lower than in the LPSF group (Fig. 7c) (log in a similar protocol used to establish the single cell fold change − 1.0, t-test p value 0.0003). These data clones [30]. Data from SCC25 were also included as a suggest that EGFR overexpression is associated with bet- reference. It is important to note that the treatment time ter response to cetuximab, consistent with the mechan- for the SCC1 and 1CC8 pair is on the order of hours vs ism of action of the therapeutic. The opposite was weeks, as used to generate the time course data. By pro- observed for FGFR1, with overexpression in SPFS vs jecting these data into the CoGAPS signatures, the rela- LPSF (Fig. 7d, log fold change 0.9, t-test p value 0.003). tionship between the sensitive SCC1 and resistant 1CC8 However, the Bossi et al. study [42] lacks DNA methyla- recapitulates the relationship between PBS and CTX tion data to assess whether FGFR1 was epigenetically reg- time course generations, respectively, in treatment ulated in these samples. Most patients with SPFS in this driven signatures (Fig. 8a, b, e, f). Conversely, CoGAPS dataset also had intrinsic resistance to cetuximab, instead signatures related to culture specific conditions failed to of acquired resistance studied in our in vitro model. produce meaningful differences between the lines (Fig. FGFR1 expression (log2 RSEM) Relative expression ( Ct) Stein-O’Brien et al. Genome Medicine (2018) 10:37 Page 14 of 22 FGFR1 and EGFR gene expression in HNSCC FGFR1 gene expression and DNA methylation ab and normal samples (TCGA) in HNSCC and normal samples (TCGA) Pearson r=-0.06 Pearson r=-0.32 p=0.33 p<0.0001 10 12 14 16 0.0 0.2 0.4 0.6 0.8 1.0 Normal Normal EGFR expression (log2 RSEM) FGFR1 methylation (beta - cg20148210) HNSCC HNSCC EGFR gene expression in HNSCC with FGFR1 gene expression in HNSCC with c d LPFS and SPFS LPFS and SPFS LPFS SPFS LPFS SPFS (n=14) (n=26) (n=14) (n=26) Fig. 7 FGFR1 gene expression and DNA methylation patterns were confirmed in independent HNSCC tumor sample datasets. a Scatter plot of gene expression for EGFR and FGFR1 in HPV-negative HNSCC samples from TCGA demonstrated that only a few HNSCC cases present increased levels of both genes and that there is no significant correlation between the expression of both genes concomitantly. b DNA methylation of FGFR1 was anti-correlated with expression in HPV-negative HNSCC TCGA samples, suggesting that upregulation of FGFR1 might be a result of promoter demethylation in primary tumors. c EGFR expression was significantly overexpressed in a group of HNSCC patients with LPFS relative to patients with SPFS in gene expression data from Bossi et al. d FGFR1 was significantly overexpressed in patients with SPFS relative to patients with LPFS in this same dataset. LPFS - long-progression-free survival; SPFS - short-progression-free survival 8c, d). Projections of the expression patterns from were measured for a prolonged exposure (11 weeks) to a CoGAPS into the cell line gene expression data were targeted therapeutic agent. Using our novel robust time also anti-correlated with projections of the methylation course integrated analysis approach, we characterized signatures in these same data. Gene expression data the molecular alterations during the development of ac- from HNSCC tumors from patients before their treat- quired cetuximab resistance using a HNSCC in vitro ment with cetuximab described in Bossi et al. [42] were model. Cell proliferation, gene expression, and DNA also analyzed. Projection into both the CoGAPS signa- methylation high-throughput analysis were performed tures of resistance and therapeutic response significantly weekly in equivalent cultures (cetuximab and PBS con- −5 stratified LPFS vs SPFS (p value = 5.2 × 10 and 3.1 × trol generations) as resistance developed. Over the −3 10 , respectively, (Fig. 8g, h). Conversely, projection in course of 11 weeks, it was possible to compare treated to the CoGAPS signature associated with culturing was (CTX) and untreated (PBS) cells grown under the same not significant (p value = 0.50, Fig. 8i). conditions. Applying robust bioinformatics algorithms, we discriminated changes associated with acquired re- Discussion sistance from those related to adaptive response to the Analysis of course omics data enables separation of cell culturing process and treatment. The SCC25 cell immediate treatment response and technical artifacts line model was chosen since this is one of the only two from molecular mechanisms of acquired therapeutic HNSCC cell lines previously used to generate isogenic resistance cetuximab-resistant cell lines [10]. However, this is the Although numerous short time course genomics studies first study to our knowledge to enable characterization of therapeutic response have been performed [49–51], of the transcriptional and epigenetic dynamics at the this is the first time that genetic and epigenetic changes early phases of therapeutic resistance. These phases EGFR expression (ILMN_1696521) FGFR1 expression (log2 RSEM) 678 9 10 11 12 6 7 8 9 10 11 12 13 FGFR1 expression (ILMN_2398261) FGFR1 expression (log2 RSEM) 6789 10 11 12 8 9 10 11 12 Stein-O’Brien et al. Genome Medicine (2018) 10:37 Page 15 of 22 HNSCCCellline gene expression projected into CoGAPS expression signatures 0.080 a b c 0.070 0.090 0.075 0.068 Response Response Response Res Res Res 0.085 Sen Sen Sen 0.070 Treatment Treatment Treatment 0.066 0.080 0 0 0 1 1 1 P-value P-value P-value * 2.1e-04 * 5.7e-02 * .93 0.065 0.064 0.075 # 1.5e-05 # 4.1e-03 # .83 SCC25 SCC1 1CC8 SCC25 SCC1 1CC8 SCC25 SCC1 1CC8 CellLine CellLine CellLine HNSCCCellline gene expression projected into CoGAPS metylation signatures 24.6 d e f * −4.4 24.4 −0.50 Response Response Response −4.5 Res 24.2 Res Res Sen Sen Sen −0.54 Treatment OG2 Treatment −4.6 Treatment 24.0 0 0 0 1 1 1 −0.58 P-value 23.8 P-value −4.7 P-value * .47 * 7.8e-03 * 1.9e-03 # .56 # 1.3e-03 # 1.6e-03 23.6 −4.8 SCC25 SCC1 1CC8 SCC25 SCC1 1CC8 SCC25 SCC1 1CC8 CellLine CellLine CellLine Gene expression from relapsed HNSCC patients from Bossi, et al projected into CoGAPS expression signatures g h i 0.14 0.130 0.17 0.125 p−value = 0.500447 0.13 p−value = 0.003070 0.16 p−value = 0.000052 0.120 0.15 0.12 0.115 0.14 0.11 0.110 0.13 0.10 Long Short Long Short Long Short Progression Free Survival length Progression Free Survival length Progression Free Survival length Fig. 8 CoGAPS gene signatures confirmed in independent in vitro and in vivo expression data. a–c Box plots of sample weights for HNSCC cetuximab sensitive cell lines SCC25, SCC1, and resistant 1CC8 in CoGAPS expression signatures. d–f Box plots of sample weights for HNSCC cetuximab-sensitive cell lines SCC25, SCC1, and resistant 1CC8 in CoGAPS methylation signatures. g-i Box plots of tumor gene expression profiles from relapsed HNSCC patients projected into CoGAPS expression signatures. Patient samples are stratified by length of progression-free survival Pattern.1 Pattern.1 Projected Pattern 1 Pattern.2 Pattern.2 Projected Pattern 2 Pattern.3 Pattern.3 Projected Pattern 3 Stein-O’Brien et al. Genome Medicine (2018) 10:37 Page 16 of 22 cannot be measured in patients due to the complexity of albeit to a lesser degree than lowly expressed genes. Be- early detection of resistance and obtaining repeated bi- cause the PatternMarker statistics includes genes that opsy samples. are uniquely associated with each inferred biological process, these highly expressed genes would be filtered Unsupervised time course analysis with CoGAPS quantifies, from associations with the dynamic conditions. To in- visualizes, and enables functional analysis of the dynamics clude these genes in the signatures defined in this study, of acquired therapeutic resistance across omics data EP5 was filtered from the calculation of the PatternMar- modalities ker statistic. Such filtering process is not required for Determining the dynamics of the molecular alterations heatmap-based visualization and filtering of flat patterns responsible for resistance requires integrated, time is recommend when defining gene signatures containing course bioinformatics analysis to quantify these alter- sets of genes that are most strongly associated with dy- ations. Based upon previous performance of Bayesian, namic changes. Patterns that reflect technical artifacts in non-negative matrix factorization algorithms in inferring the data, such as EP4, should be retained in the Pattern- dynamic regulatory networks for targeted therapeutics Marker analysis to limit the signatures associated with [49, 50], we selected CoGAPS [23] for analysis of gene inferred processes to retain only biologically relevant expression data from our time course experiment. genes. We note that these PatternMarker statistic are CoGAPS have already proven highly effective in relating similar to the D-scores proposed in Zhu et al. [52] and gene expression changes to patterns related to EGFR in- that application of this statistic may require similar fil- hibition [41], perturbation of nodes in the EGFR net- tering to retain highly expressed genes. In the case of work [39], and time course dynamics of targeted DNA methylation, these PatternMarker genes also in- therapeutics. In this dataset, CoGAPS analysis of gene clude genes representing driver alterations in resistance. expression data from cetuximab-resistant clones distin- The DNA methylation data did not require filtering guished the patterns for immediate gene expression when applying the PatternMarker statistic since no flat changes and patterns for long-term changes associated pattern was detected. However, transcriptional regula- with acquired resistance. Gene expression signatures for tion by epigenetic alterations or in pathways involves resistance to EGFR inhibitors in two additional cell lines simultaneous co-regulation of multiple genes. This co- (one HNSCC and one NSCLC) from previous studies regulation is reflected in the reuse of genes in CoGAPS [38, 41] were significantly enriched in both types of gene weights associated with each pattern. Therefore, es- CoGAPS patterns. Since these previous resistance signa- timates of pathway dynamics from transcriptional data tures were learned from case-control studies without require accounting for all genes with gene set enrich- multiple time point measurements, we concluded that ment statistics instead of the PatternMarker statistic. our time course data are instrumental in discriminating Thus, we hypothesize that the PatternMarker statistic is the signatures of immediate therapeutic response from robust for visualization and biomarker identification. On signatures of acquired resistance. the other hand, gene set enrichment of the CoGAPS In spite of the complexities of the data integration, the gene weights corresponding to each pattern and stored weight of each sample in patterns inferred by CoGAPS in the amplitude matrix are essential for characterization reflects the dynamics of the process in each data modal- of functional alterations in pathways. ity. These patterns are learned completely unsupervised from the data and do not require any gene selection or Integrated genomics analysis of time course data comparison between time points relative to any refer- demonstrates that DNA methylation changes follow ence control. The CoGAPS analysis of the time course transcriptional changes, leading to the model where data demonstrates that applying matrix factorization al- methylation stabilizes the resistance phenotype gorithms for genomics can reconstruct signals associated Collecting treated and untreated cells to obtain paired with phenotypes from time course, omics data. The measurements of methylation and gene expression en- genes associated with CoGAPS patterns had weights that abled us to evaluate whether changes in DNA methyla- were non-zero in multiple patterns. The PatternMarker tion impact gene expression. Including a PBS control at statistic [26] enable further selection of the genes that every time point also enabled the discrimination of the are uniquely associated with each pattern. Creating a changes that result from an adaptive response to therapy heatmap of the genomics profiles for these genes en- from changes that result from maintaining cells in cul- abled novel, heatmap-based visualization of the temporal ture. CoGAPS analysis of DNA methylation data denotes dynamics in the omics data. CoGAPS analysis of gene only changes associated with acquired resistance, in con- expression data contains a flat pattern (EP5), which in- trast to the immediate expression changes observed with cludes all highly expressed genes. These genes may also cetuximab treatment. Thus, while therapeutic response change in association with the experimental conditions, can drive massive changes in gene expression, only the Stein-O’Brien et al. Genome Medicine (2018) 10:37 Page 17 of 22 subset of expression changes associated with the devel- gene regulation, associated with acquired therapeutic resist- opment of resistance have corresponding epigenetic sig- ance. Future investigation into the chromatin remodeling natures, suggesting that the methylation landscape is mechanisms will also test whether chromatin alterations important for the development of acquired resistance. follow the changes in expression and occur in combination The CoGAPS patterns in gene expression that are asso- with altered methylation patterns to drive epigenetic regula- ciated with acquired cetuximab resistance gradually tion of resistance. change over the time course (EP1 and EP2). On the other hand, the CoGAPS patterns for DNA methylation changes have a sharp transition at the generation at Time course data are essential to distinguish clonal which resistance is acquired (CTX-G4). These patterns outgrowth from transcriptional rewiring that give rise to (MP2 and MP3) reflect a delayed but more rapid change stable acquired resistance to therapy in DNA methylation. The time delays between alter- Besides the immediate changes in gene expression ations in DNA methylation and gene expression pose a followed by the gradual methylation switch, it is also in- further computational challenge for integrated, time teresting to note these effects in the proliferation rates course genomics analyses. The vast majority of inte- during the 11 weeks of treatment. Initially, the prolifera- grated analysis algorithms assume one-to-one mapping tion of the population of cetuximab-sensitive cells is of genes in different data platforms or seek common pat- slower when compared to the untreated controls, reflect- terns or latent variables across them [53]. Such ap- ing therapy effectiveness. Early and progressively, the proaches would fail to capture the early changes from cells develop molecular changes to overcome the EGFR cetuximab treatment that impact only gene expression; blockade. However, this process starts in just a small time delays between DNA methylation and gene expres- number of clones and the increase in proliferation is still sion patterns and different gene usage in each pattern. It not enough to surpass the growth rate of the untreated is essential to develop new integrated algorithms to sim- cells. As soon as the population of resistant cells is larger ultaneously distinguish both patterns that are shared than the number of sensitive cells, the proliferation rate across data types and that are unique to each platform. is now higher than in the untreated controls. At some For time course data analysis, these algorithms must also point, we observe the stabilization of the proliferation model regulatory relationships that may give rise to tim- rates in the cetuximab-treated cells, probably due to the ing delays, such as epigenetic silencing of gene expres- fact the culture is now dominated by the population of sion. However, as we observed with the unanticipated resistant clones. Although stable, the proliferation rate changes in DNA methylation following and not preced- of these resistant clones is significantly higher than that ing gene expression, they must also consider delays of the untreated cells. This increased proliferation rate is resulting from larger phenotypic changes such as the consistent with the rapid increase in tumor volume ob- stability of the therapeutic resistant phenotype. served clinically once patients develop resistance to ther- The relative timing of change in DNA methylation and apy. Tracing the increase in the population of resistant gene expression is consistent with previous observations cells and their proliferation rates in vivo requires novel that gene expression changes precede DNA methylation al- techniques to biopsy or image tumors at intermediate INK4A terations in genes critical for cancer progression. P16 time points of treatment. and GSTP1 are tumor suppressor genes for which In a recent study, gene expression changes were asso- transcription silencing was found to occur prior to DNA ciated with a transient resistant phenotype present in hypermethylation and chromatin changes. The temporal melanoma cell lines before vemurafenib administration delay observed between expression and methylation [56]. Once the melanoma cells were exposed to the drug, patterns in our time course provides transcriptome-wide additional changes in gene expression are detected and evidence of this phenomenon. Specifically, that epigenetic are later accompanied by changes in chromatin structure changes are necessary to stabilize gene expression aberrant [56]. These findings, together with our time course ob- profile and will be followed by modification into a silenced servations, suggest that in the heterogeneous tumor en- methylation state, resulting in tumor progression [54, 55]. vironment the existence of some cells expressing specific Our integrated RNA-seq and DNA methylation analysis marker genes can trigger cellular reprogramming as corroborates the fact that gene expression changes occur soon as the targeted therapy is initiated. Upon drug ad- earlier to epigenetic alterations and suggests that DNA ministration, the number of genes with aberrant expres- methylation is essential to maintain the changes in gene ex- sion increases and is followed by other epigenetic and pression in this acquired cetuximab resistance model. genetic changes that will shift the transient resistant Additional time course data tracing other in vitro and in state into a stable phenotype. This finding on acquired vivo models of HNSCC are essential to generalize the rela- resistance development could dramatically change the tive timing of molecular changes, and thus mechanisms of course of treatment with targeted therapeutic agents. Stein-O’Brien et al. Genome Medicine (2018) 10:37 Page 18 of 22 The precise characterization of resistant gene signatures inhibition resulted in phospho-ERK decreased expres- and their timing are crucial to determine the correct sion that was restored when FGF2 was added to the point during the patients’ clinical evolution to introduce culture, suggesting that an autocrine FGF-FGFR path- alternative therapeutic strategies. This way, secondary way is one of the mechanisms of resistance to gefitinib. interventions would start before the stable resistant However, the cell lines evaluated in both studies were phenotype is spread among the tumor cells resulting in intrinsically resistant to gefitinib. In our model, we in- prolonged disease control and substantial increase in duced resistance to cetuximab and observed FGFR1 overall survival. gain of expression and significant anti-correlation with the DNA methylation. We additionally evaluated the expression of other FGFRs and FGFs that were identi- DNA methylation changes in FGFR1 are associated with fied by the PatternMarker statistic. Although it is not changes in FGFR signaling found as a PatternMarker in our analysis pipeline, Among the genes we identified with the canonical rela- FGF2 is upregulated in the cetuximab generations when tionship between expression and methylation, FGFR1 compared to the PBS generations as observed with present with increased gene expression accompanied by FGFR1 (Additional file 1: Figure S14). Thus, our data loss of CpG methylation. FGFR1 is a receptor tyrosine corroboratesand extendsthispreviousevidencefrom kinase that regulates downstream pathways, such as intrinsically resistant lines that one of the mechanisms PI3K/AKT, and RAS/MAPK, which are also regulated driving resistance to EGFR inhibition is the FGF-FGFR by EGFR [57]. Its overexpression has been previously autocrine pathway. This observation adds another evi- associated with resistance to EGFR inhibitors in other dence that the computational approach used in this cancer types including HNSCC [27–29]. To our know- study is robust once it is capable of identifying mecha- ledge this is the first study showing that epigenetic al- nisms previously described in other models resistant to terations are associated with changes in FGFR1 EGFR blockade. expression in HNSCC during the development acquired Our previously developed bioinformatics algorithms cetuximab resistance. FGFR1 upregulation combined for the identification of gene expression and epigenetic with promoter hypomethylation was previously de- patterns progression over time proved to be consistent, scribed in rhabdomyosarcomas [58]. Other studies de- since they also detected canonical changes found to be scribed that FGFR1 increased levels is a common driving this mechanism among innumerous new poten- feature in different tumor types, such as glioblastoma tial candidates for acquired resistance. The integrated [59] and cancers of the breast [60], lung [61], prostate computational analysis was possible due to an experi- [62], bladder [63], ovarian [64], colorectal [27], and mental approach developed to account for molecular HNSCC [29, 65, 66]. FGFR1 is involved in resistance changes due to adaptive responses to the culturing sys- mechanisms against EGFR inhibitors [27, 46–48], such tem and the immediate addition of cetuximab. Here, we as cetuximab and gefitinib. Together, the TCGA and present a novel integrated analysis protocol to evaluate Bossi et al. dataset analyses corroborate our findings molecular changes measured by different high through- that FGFR1 gene expression is regulated by epigenetic put techniques over a prolonged time of treatment with changes in HNSCC. Altogether, the epigenetic alter- an FDA-approved targeted therapeutic agent. The lack ation of FGFR1 represents a candidate biomarker of re- of in vivo experimentation to validate our findings was sistance to cetuximab and further studies are critical to compensated by the analysis of two public datasets of identify combination therapies for HNSCC patients that HNSCC, showing that our in vitro findings were also develop acquired cetuximab resistance. present in patients’ samples. Our findings, together with The increased levels of FGFRs and FGFs are believed Marshall et al. [47] and Marek et al. [46], are a strong to play a role in an autocrine mechanism in HNSCC evidence that FGFR1 plays a crucial role in a significant and NSCLC cell lines with intrinsic resistance to the proportion of cases that are resistant to cetuximab or EGFR inhibitor, gefitinib. Using publicly available gene gefitinib. The translational implications are notable expression microarray datasets, Marshall et al. [47]and since FGFR1 inhibition can be used in combination Marek et al. [46] verified concomitant increased levels with EGFR blockade to retard acquired resistance or of FGFRs and their specific FGFs ligands. Particularly, overcome intrinsic resistance. It is important to men- FGFR1 and FGF2 upregulation was observed in the tion that FGFR inhibitors are being currently evaluated same resistant cell lines and hypothesized to be the by clinical trials and could soon become a potential mechanism behind resistance. This was corroborated new therapeutic option for many cancer patients [57]. by functional experiments showing that cells treated Future work evaluating how these combinations impact with pan-FGFR inhibitor were less prone to anchorage- thetimingofacquiredresistance are essential to independent growth. Also, FGF2 silencing or FGFR1 determine the molecular mechanisms that shift Stein-O’Brien et al. Genome Medicine (2018) 10:37 Page 19 of 22 dominant signaling pathways in cancer and thereby when treated with cisplatin [68]. These data and the in- drive resistance. trinsic sensitivity of SCC25 to EGFR inhibition suggests that therapeutic resistance results from random selection of a pre-existing resistant clone. The heterogeneity in CoGAPS gene signatures from the SCC25 time course methylation profiles reflected the complexity of the resist- model are associated with molecular changes in ance mechanisms that can arise from combination therap- additional in vitro models and human tumor data ies in heterogeneous tumors. Future work extending these The main limitation of the current study was the use of protocols to in vivo models is essential to determine the a single cell line model. SCC25 is intrinsically sensitive role of the microenvironment in inducing therapeutic to cetuximab and from this single cell line model, we resistance. Developing in vivo models with acquired thera- generated two groups of samples (CTX and PBS genera- peutic resistance presents numerous technical challenges tions) over the course of 11 weeks. High-throughput that must first be addressed before such time course measurements and analysis were performed for a total of sampling is possible [9]. Pinpointing precise molecular 22 samples. The collection of multiple data points in the predictors of therapeutic resistance will facilitate the analysis had to be accounted for when determining the identification of unprecedented biomarkers and reveal the number of cell lines to be included in the study. We mechanisms by which to overcome acquired therapeutic nonetheless compared our data to gene signatures from resistance to most therapies used to treat cancer. the other isogenic HNSCC resistant model 1CC8 [10], an independent resistant model to an EGFR inhibitor in NSCLC [38], and human tumor data from HNSCC pa- Conclusions tients before cetuximab treatment [42]. Besides the By developing a novel bioinformatics pipeline for inte- number of samples, we also had to take into consider- grated time course analysis, we measured the changes ation the potential batch and technical effects of broad in gene expression and DNA methylation during the cross-platform profiling. Nevertheless, the analysis of progression from an intrinsic cetuximab responsive pretreatment HNSCC patient samples from TCGA [36] state to the acquired resistant phenotype using an in and another study [42] confirmed that our finding that vitro HNSCC cell line model. Specifically, this pipeline FGFR1 is upregulated and demethylated in HNSCC and includes: (1) CoGAPS analysis of each platform inde- associated with acquired resistance to cetuximab is also pendently; (2) gene selection with the PatternMarker a mechanism involved in intrinsic resistance to the tar- statistic for visualization and CoGAPS gene set analysis geted therapy. of the CoGAPS gene profiles for pathway analysis; (3) comparisons of patterns to known phenotypes infer Transcriptional heterogeneity is critical in acquired their relative timing; (4) anti-correlation between DNA therapeutic resistance, and future time course single cell methylation patterns and gene expression to infer epi- data will pinpoint precise molecular predictors of genetically regulated genes; and (5) evaluation of Pat- therapeutic resistance ternMarker genes and projection of the CoGAPS gene The in vitro protocol for time course sampling developed profiles to learn relevance of inferred gene signatures in in this study has the additional advantage of aggregating new datasets. This pipeline revealed massive changes in potentially heterogeneous mechanisms of resistance in- gene expression and identified and discriminated the creasing the signal of changes in any cetuximab-resistant different patterns associated with resistance or cell cul- sub-clone. For example, we observed demethylation and turing conditions. This analysis demonstrates that com- overexpression of FGFR1 in the pooled cells, but only a pressed sensing matrix factorization algorithms can single stable clone generated from the same SCC25 cell identify gene signatures associated with the dynamics line in a previous study (CTXR10) had upregulation of of phenotypic changes from genomics data collected FGFR1 [31]. This finding suggests that tumor heterogen- over the time course. In this case, the gene expression eity also plays a role in acquired resistance to targeted patterns relevant to resistance were later followed by therapies and enables different pathways to be used to epigenetic alterations. Our main conclusion is that bypass the silenced target within the same tumor. using our bioinformatics approach we are able to deter- Heterogeneity of SCC25 cetuximab-resistant clones has mine that the resistant phenotype is driven by gene ex- been observed previously [31]. Recent single cell RNA-seq pression changes that would confer the cancer cells data of SCC25 has shown that there is considerable tran- adaptive advantages to the treatment with cetuximab. scriptional heterogeneity in this cell line before treatment Finally, the integrated analysis show that the stability of [67]. Other cancer therapies are influenced by heterogen- the resistant state is dependent on epigenetic changes eity and outgrowth of resistant clones, as was observed in that will make these new gene signatures heritable to single cell clones isolated from the HNSCC cell line FaDu expand the phenotype to the daughter cells. The Stein-O’Brien et al. Genome Medicine (2018) 10:37 Page 20 of 22 bioinformatics pipeline we developed is also significant acid; RNA-seq: Ribonucleic acid sequencing; rRNA: Ribosomal ribonucleic acid; SPFS: Short-progression-free survival; STR: Short tandem repeat; to clinical practice, since it pointed the time course of TCGA: The Cancer Genome Atlas; TFAP2A: Transcription factor AP-2 alpha molecular changes associated with acquired cetuximab resistance and suggests that the resistant phenotype Acknowledgements can be reversed if alternative interventions are intro- We thank JHMI Deep Sequencing & Microarray Core and SKCCC Microarray duced before epigenetic alterations to the genes driving Core Facility on performing and providing advice on RNA-seq and DNA methylation hybridization arrays, respectively; S. Boca, B. Kerr, S. Floor, C. Mak, acquired resistance. Moreover, the computational T. Ou, D. Sidransky, L. M. Weiner, F. Zamuner, K. Zambo, and members of approach we describe here can be applied to time NewPISlack for critical comments and feedback during the preparation of course studies using other tumor type models and the manuscript. In memory of Elizabeth Klein. targeted therapies. Funding This work was supported by NIH Grants R01CA177669, R21DE025398, Additional files P30CA006973, R01DE017982, SPORE P50DE019032, and the Johns Hopkins University Catalyst Award. Additional file 1: Supplemental Figures S1 - S14. Figure S1 - Time course approach to induce resistance to cetuximab and measure gene Availability of data and materials expression and DNA methylation changes; Figure S2 - Anchorage- Unless otherwise specified, all genomics analyses were performed in R and independent growth of cetuximab generation 10 (CTX-G10); Figure S3 - code for these analyses is available from https://sourceforge.net/projects/ Heatmap and hierarchical clustering of gene expression values in 11 gen- scc25timecourse. All transcriptional and epigenetic data of the cell lines from erations of SCC25 cells treated with PBS as control (black columns) and this paper have been deposited in GEO (GSE98815). DNA methylation data with 100 nM of cetuximab (red columns) to acquire resistance.; Figure of head and neck tumors after treatment with cetuximab are available in S4 - Time course gene expression comparison to previously known gene GEO (GSE110996). signatures of resistance to EGFR inhibitors; Figure S5 - Expected gene ex- pression values for genes in each CoGAPS pattern inferred from gene ex- pression data over generations of PBS control (black lines) or treatment Authors’ contributions with 100 nM of cetuximab (red lines); Figure S6 - Heatmap of gene set GS, LTK, SL, CHC, and EJF planned, designed and wrote the manuscript with analysis scores for targets of transcription factors in theEGFR network, tar- input from all authors. GS, LTK, SL, MT, CHC, and EJF contributed to the gets of the AP-2alpha transcription factors associated with cetuximabre- development of methodology. SS provided DNA from HNSCC patients sponse, and cetuximab resistance signatures in CoGAPS patterns; Figure following cetuximab treatment for DNA methylation analysis. GS, SL, and EJF S7 - Heatmaps of Pearson correlation coefficients between CoGAPS gene performed analysis and interpretation of data (e.g. computational analysis). expression and DNA methylation patterns; Figure S8 - Gene expression RR, HO, HC, MC, AF, LVD, SS, JA, and DAG participated in development of heatmap for the time course experiment vs. single cell resistant clones methodology and provided technical and material support. RR, LVD, EI, and experiment; Figure S9 - DNA methylation heatmap for the time course DAG participated in the review and/or revision of the manuscript. All authors experiment vs. single cell resistant clones experiment; Figure S10 - Mi- discussed the data and contributed to the manuscript preparation. CHC and croscopy images of cetuximab single cell resistant clones (CTXR); Figure EJF instigated and supervised the project. All authors read and approved the S11 - Epigenetically regulated pattern marker genes associated with re- final manuscript. sistance presenting significant anti-correlation between gene expression and DNA methylation in thecetuximab single cell resistant clones (CTXR); Figure S12 - Proliferation assay of cetuximab resistant single cell clones Ethics approval and consent to participate (CTXR); Figure S13 - Bar plot of FGFR1 expression in human tumor sam- HNSCC samples and patient information collection were approved by the ples pre (black) and post (red) treatment with cetuximab; Figure S14 Independent Ethics Committee and the Belgian Health Authorities (ID: 2008/ - Heatmap of FGF family genes methylation (A) and expression (B) in time 20MAR/090) and were conducted in accordance with the Declaration of coursedata for treated and control samples. (PDF 22593 kb) Helsinki (October 2000). It was prospectively planned to perform translational research and patients gave their informed consent for repeated biopsies. Additional file 2: Table S1. PatternMarker genes for the expression and methylation CoGAPS patterns and list of genes with significant anti- correlation between gene expression and methylation. (XLSX 117 kb) Competing interests Additional file 3: Table S2. P values of MSigDB hallmark pathways for The authors declare that they have no competing interests. CoGAPS expression patterns calculated with the permutation-based statistic in calcCoGAPSStat. (XLSX 99 kb) Additional file 4: Table S3. P values of MSigDB hallmark pathways Publisher’sNote for CoGAPS DNA methylation patterns calculated with the permutation- Springer Nature remains neutral with regard to jurisdictional claims in published based statistic in calcCoGAPSStat. (XLSX 79 kb) maps and institutional affiliations. Author details Abbreviations Institute of Genetic Medicine, Johns Hopkins University, Baltimore, MD, USA. AKT: AKT serine/threonine kinase; ATCC: American Type Culture Collection; Department of Oncology, Sidney Kimmel Comprehensive Cancer Center, CoGAPS: Coordinated Gene Activity in Pattern Sets; CTX: Cetuximab; Johns Hopkins University, Baltimore, MD, USA. Department of Head and CTXR: Single cell cetuximab resistant clone; DNA: Deoxyribonucleic acid; Neck-Endocrine Oncology, Moffitt Cancer Center, Tampa, FL, USA. EGFR: Endogenous growth factor receptor; FDA: Food and Drug Department of Otorhinolaryngology-Head and Neck Surgery, Keio University Administration; FGFR1: Fibroblast growth factor receptor 1; GEO: Gene School of Medicine, Tokyo, Japan. Department of Surgery - Otolaryngology– Expression Omnibus; GSTP1: Glutathione S-Transferase pi 1; Head and Neck Surgery, University of Utah, |Salt Lake City, UT, USA. Head GWCoGAPS: Genome Wide Coordinated Gene Activity in Pattern Sets; and Neck Surgery Unit, St Luc University Hospital, Brussels, Belgium. HNSCC: Head and neck squamous cell carcinoma; HRAS: HRAS proto- Laboratory of Systems Biology and Computational Genetics, Vavilov Institute oncogene; JHMI: Johns Hopkins Medical Institutions; LPFS: Long-progression- of General Genetics, Russian Academy of Sciences, Moscow, Russia. free survival; MAPK: Mitogen activated kinase-like protein; NSCLC: Non-small- Department of Surgery, UC San Diego Moores Cancer Center, La Jolla, CA, cell lung cancer; PBS: Phosphate buffered saline; PI3K: Phosphoinositide-3- USA. Department of Otolaryngology-Head and Neck Surgery, Johns Hopkins kinase regulatory subunit 1; RIN: RNA integrity number; RNA: Ribonucleic University, Baltimore, MD, USA. Stein-O’Brien et al. Genome Medicine (2018) 10:37 Page 21 of 22 Received: 23 January 2018 Accepted: 1 May 2018 24. Ochs MF, Fertig EJ. Matrix factorization for transcriptional regulatory network inference. IEEE Symp Comput Intell Bioinforma Comput Biol Proc. 2012;2012:387–96. 25. Fertig EJ, Markovic A, Danilova LV, Gaykalova DA, Cope L, Chung CH, et al. Preferential activation of the hedgehog pathwaybyepigenetic modulations in References HPV negative HNSCC identified with meta-pathway analysis. PLoS One. 2013; 1. Sawyers C. Targeted cancer therapy. Nature. 2004;432:294–7. e78127:8. 2. Hyman DM, Taylor BS, Baselga J. Implementing Genome-Driven Oncology. 26. Stein-O’Brien GL, Carey JL, Lee WS, Considine M, Favorov AV, Flam E, et al. Cell. 2017;168:584–99. PatternMarkers & GWCoGAPS for novel data-driven biomarkers via whole 3. Engelman JA, Zejnullahu K, Mitsudomi T, Song Y, Hyland C, Park JO, et al. transcriptome NMF. Bioinformatics. 2017;33:1892–4. MET amplification leads to gefitinib resistance in lung cancer by activating 27. Azuma K, Kawahara A, Sonoda K, Nakashima K, Tashiro K, Watari K, et al. FGFR1 ERBB3 signaling. Science. 2007;316:1039–43. activation is an escape mechanism in human lung cancer cells resistant to 4. Hata AN, Niederst MJ, Archibald HL, Gomez-Caraballo M, Siddiqui FM, afatinib, a pan-EGFR family kinase inhibitor. Oncotarget. 2014;5:5908–19. Mulvey HE, et al. Tumor cells can follow distinct evolutionary paths to 28. Bertotti A, Papp E, Jones S, Adleff V, Anagnostou V, Lupo B, et al. The become resistant to epidermal growth factor receptor inhibition. Nat Med. genomic landscape of response to EGFR blockade in colorectal cancer. 2016;22:262–9. Nature. 2015;526:263–7. 5. Misale S, Yaeger R, Hobor S, Scala E, Janakiraman M, Liska D, et al. 29. Koole K, Brunen D, van Kempen PMW, Noorlag R, de Bree R, Emergence of KRAS mutations and acquired resistance to anti-EGFR therapy Lieftink C, et al. FGFR1 is a potential prognostic biomarker and therapeutic in colorectal cancer. Nature. 2012;486:532–6. target in head and neck squamous cell carcinoma. Clin Cancer Res. 6. Pietrantonio F, Vernieri C, Siravegna G, Mennitto A, Berenato R, Perrone F, et al. 2016;22:3884–93. Heterogeneity of acquired resistance to anti-EGFR monoclonal antibodies in 30. Hatakeyama H, Cheng H, Wirth P, Counsell A, Marcrom SR, Wood CB, et al. patients with metastatic colorectal cancer. Clin Cancer Res. 2017;23:2414–22. Regulation of heparin-binding EGF-like growth factor by miR-212 and 7. Vincenzi B, Zoccoli A, Pantano F, Venditti O, Galluzzo S. Cetuximab: from acquired cetuximab-resistance in head and neck squamous cell carcinoma. bench to bedside. Curr Cancer Drug Targets. 2010;10:80–95. PLoS One. 2010;e12702:5. 8. Boeckx C, Weyn C, Vanden Bempt I, Deschoolmeester V, Wouters A, 31. Cheng H, Fertig EJ, Ozawa H, Hatakeyama H, Howard JD, Perez J, et al. Specenier P, et al. Mutation analysis of genes in the EGFR pathway in Head Decreased SMAD4 expression is associated with induction of epithelial-to- and Neck cancer patients: implications for anti-EGFR treatment response. mesenchymal transition and cetuximab resistance in head and neck BMC Res Notes. 2014;7:337. squamous cell carcinoma. Cancer Biol Ther. 2015;16:1252–8. 9. Quesnelle KM, Wheeler SE, Ratay MK, Grandis JR. Preclinical modeling of 32. Wang K, Singh D, Zeng Z, Coleman SJ, Huang Y, Savich GL, et al. MapSplice: EGFR inhibitor resistance in head and neck cancer. Cancer Biol Ther. 2012; Accurate mapping of RNA-seq reads for splice junction discovery. Nucleic 13:935–45. Acids Res. 2010;38:e178. 10. Wheeler DL, Huang S, Kruser TJ, Nechrebecki MM, Armstrong EA, Benavente 33. Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data S, et al. Mechanisms of acquired resistance to cetuximab: role of HER (ErbB) with or without a reference genome. BMC Bioinformatics. 2011;12:323. family members. Oncogene. 2008;27:3944–56. 34. Fortin J-P, Labbe A, Lemire M, Zanke BW, Hudson TJ, Fertig EJ, et al. 11. Narayan M, Wilken JA, Harris LN, Baron AT, Kimbler KD, Maihle NJ. Functional normalization of 450k methylation array data improves Trastuzumab-induced HER reprogramming in “resistant” breast carcinoma replication in large cancer studies. Genome Biol. 2014;15:503. cells. Cancer Res. 2009;69:2191–4. 35. Bidaut G. Interpreting and comparing clustering experiments through graph 12. Smyth GK. Linear models and empirical Bayes methods for assessing visualization and ontology statistical enrichment with the ClutrFree package. differential expression in microarray experiments. Stat Appl Genet Mol Biol. In: Ochs MF, Casagrande JT, Davuluri RV, editors. Biomed. Inform. Cancer 2004;3:Article3. Res. Boston, MA: Springer US; 2010. p. 315–33. http://link.springer.com/10. 13. Ernst J, Nau GJ, Bar-Joseph Z. Clustering short time series gene expression 1007/978-1-4419-5714-6_19. Accessed 23 Jan 2018. data. Bioinforma Oxf Engl. 2005;21(Suppl 1):i159–68. 14. Lin D, Shkedy Z, Yekutieli D, Burzykowski T, Göhlmann HWH, De Bondt A, et 36. Lawrence MS, Sougnez C, Lichtenstein L, Cibulskis K, Lander E, Gabriel SB, al. Testing for trends in dose-response microarray experiments: a et al. Comprehensive genomic characterization of head and neck squamous comparison of several testing procedures, multiplicity and resampling-based cell carcinomas. Nature. 2015;517:576–82. inference. Stat Appl Genet Mol Biol. 2007;6:Article26. 37. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: a knowledge-based approach for 15. Aryee MJ, Gutiérrez-Pabello JA, Kramnik I, Maiti T, Quackenbush J. An interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. improved empirical bayes approach to estimating differential gene 2005;102:15545–50. expression in microarray time-course data: BETR (Bayesian Estimation of 38. Coldren CD. Baseline gene expression predicts sensitivity to gefitinib in Temporal Regulation). BMC Bioinformatics. 2009;10:409. non-small cell lung cancer cell lines. Mol Cancer Res. 2006;4:521–8. 16. Liao JC, Boscolo R, Yang Y-L, Tran LM, Sabatti C, Roychowdhury VP. Network 39. Fertig EJ, Ren Q, Cheng H, Hatakeyama H, Dicker AP, Rodeck U, et al. Gene component analysis: reconstruction of regulatory signals in biological expression signatures modulated by epidermal growth factor receptor systems. Proc Natl Acad Sci U S A. 2003;100:15522–7. activation and their relationship to cetuximab resistance in head and neck 17. Bonneau R, Reiss DJ, Shannon P, Facciotti M, Hood L, Baliga NS, et al. The squamous cell carcinoma. BMC Genomics. 2012;13:160. Inferelator: an algorithm for learning parsimonious regulatory networks from systems-biology data sets de novo. Genome Biol. 2006;7:R36. 40. McCall MN, Bolstad BM, Irizarry RA. Frozen robust multiarray analysis (fRMA). 18. Ernst J, Vainas O, Harbison CT, Simon I, Bar-Joseph Z. Reconstructing Biostat Oxf Engl. 2010;11:242–53. dynamic regulatory maps. Mol Syst Biol. 2007;3:74. 41. Fertig EJ, Ozawa H, Thakar M, Howard JD, Kagohara LT, Krigsfeld G, et al. 19. Seok J, Xiao W, Moldawer LL, Davis RW, Covert MW. A dynamic network of CoGAPS matrix factorization algorithm identifies transcriptional changes in transcription in LPS-treated human subjects. BMC Syst Biol. 2009;3:78. AP-2alpha target genes in feedback from therapeutic inhibition of the EGFR 20. Naegle KM, Welsch RE, Yaffe MB, White FM, Lauffenburger DA. MCAM: multiple network. Oncotarget. 2016;7:73845–64. clustering analysis methodology for deriving hypotheses and insights from 42. Bossi P, Bergamini C, Siano M, Cossu Rocca M, Sponghini AP, Favales F, high-throughput proteomic datasets. PLoS Comput Biol. 2011;7:e1002119. et al. Functional genomics uncover the biology behind the responsiveness of head and neck squamous cell cancer patients to cetuximab. Clin Cancer 21. Fernández MA, Rueda C, Peddada SD. Identification of a core set of Res. 2016;22:3961–70. signature cell cycle genes whose relative order of time to peak expression is conserved across species. Nucleic Acids Res. 2012;40:2823–32. 43. Schmitz S, Bindea G, Albu RI, Mlecnik B, Machiels J-P. Cetuximab 22. Wise A, Bar-Joseph Z. SMARTS: reconstructing disease response networks promotes epithelial to mesenchymal transition and cancer associated from multiple individuals using time series gene expression data. fibroblasts in patients with head and neck cancer. Oncotarget. Bioinforma Oxf Engl. 2015;31:1250–7. 2015;6:34288–99. 23. Fertig EJ, Ding J, Favorov AV, Parmigiani G, Ochs MF. CoGAPS: an R/C++ 44. Fortin J-P, Triche TJ, Hansen KD. Preprocessing, normalization and package to identify patterns and biological process activity in transcriptomic integration of the Illumina HumanMethylationEPIC array with minfi. data. Bioinforma Oxf Engl. 2010;26:2792–3. Bioinformatics. 2017;33:558–60. Stein-O’Brien et al. Genome Medicine (2018) 10:37 Page 22 of 22 45. Mootha VK, Lindgren CM, Eriksson K-F, Subramanian A, Sihag S, Lehar J, 67. Puram SV, Tirosh I, Parikh AS, Patel AP, Yizhak K, Gillespie S, et al. Single-cell et al. PGC-1α-responsive genes involved in oxidative phosphorylation transcriptomic analysis of primary and metastatic tumor ecosystems in head are coordinately downregulated in human diabetes. Nat Genet. 2003; and neck cancer. Cell. 2017;171:1611–1624.e24. 34:267–73. 68. Niehr F, Eder T, Pilz T, Konschak R, Treue D, Klauschen F, et al. Multilayered 46. Marek L, Ware KE, Fritzsche A, Hercule P, Helton WR, Smith JE, et al. Fibroblast omics-based analysis of a head and neck cancer model of cisplatin growth factor (FGF) and FGF receptor-mediated autocrine signaling in non- resistance reveals intratumoral heterogeneity and treatment-induced clonal small-cell lung cancer cells. Mol Pharmacol. 2009;75:196–207. selection. Clin Cancer Res. 2018;24:158–68. 47. Marshall ME, Hinz TK, Kono SA, Singleton KR, Bichon B, Ware KE, et al. Fibroblast growth factor receptors are components of autocrine signaling networks in head and neck squamous cell carcinoma cells. Clin Cancer Res. 2011;17:5016–25. 48. Wynes MW, Hinz TK, Gao D, Martini M, Marek LA, Ware KE, et al. FGFR1 mRNA and protein expression, not gene copy number, predict FGFR TKI sensitivity across all lung cancer histologies. Clin Cancer Res. 2014;20:3299–309. 49. Ochs MF, Rink L, Tarn C, Mburu S, Taguchi T, Eisenberg B, et al. Detection of treatment-induced changes in signaling pathways in gastrointestinal stromal tumors using transcriptomic data. Cancer Res. 2009;69:9125–32. 50. Hafner M, Niepel M, Chung M, Sorger PK. Growth rate inhibition metrics correct for confounders in measuring sensitivity to cancer drugs. Nat Methods. 2016;13:521–7. 51. Hill SM, Heiser LM, Cokelaer T, Unger M, Nesser NK, Carlin DE, et al. Inferring causal molecular networks: empirical assessment through a community- based effort. Nat Methods. 2016;13:310–8. 52. Zhu X, Ching T, Pan X, Weissman SM, Garmire L. Detecting heterogeneity in single- cell RNA-Seq data by non-negative matrix factorization. PeerJ. 2017;5:e2888. 53. Tyekucheva S, Marchionni L, Karchin R, Parmigiani G. Integrating diverse genomic data using gene sets. Genome Biol. 2011;12:R105. 54. Bachman KE, Park BH, Rhee I, Rajagopalan H, Herman JG, Baylin SB, et al. Histone modifications and silencing prior to DNA methylation of a tumor suppressor gene. Cancer Cell. 2003;3:89–95. 55. Stirzaker C, Song JZ, Davidson B, Clark SJ. Transcriptional gene silencing promotes DNA hypermethylation through a sequential change in chromatin modifications in cancer cells. Cancer Res. 2004;64:3871–7. 56. Shaffer SM, Dunagin MC, Torborg SR, Torre EA, Emert B, Krepler C, et al. Rare cell variability and drug-induced reprogramming as a mode of cancer drug resistance. Nature. 2017;546:431–5. 57. Babina IS, Turner NC. Advances and challenges in targeting FGFR signalling in cancer. Nat Rev Cancer 2017;17:318–32. 58. Goldstein M, Meller I, Orr-Urtreger A. FGFR1 over-expression in primary rhabdomyosarcoma tumors is associated with hypomethylation of a 5’ CpG island and abnormal expression of the AKT1, NOG, and BMP4 genes. Genes Chromosomes Cancer. 2007;46:1028–38. 59. Rand V, Huang J, Stockwell T, Ferriera S, Buzko O, Levy S, et al. Sequence survey of receptor tyrosine kinases reveals mutations in glioblastomas. Proc Natl Acad Sci U S A. 2005;102:14344–9. 60. Andre F, Bachelot T, Campone M, Dalenc F, Perez-Garcia JM, Hurvitz SA, et al. Targeting FGFR with Dovitinib (TKI258): preclinical and clinical data in breast cancer. Clin Cancer Res. 2013;19:3693–702. 61. Cihoric N, Savic S, Schneider S, Ackermann I, Bichsel-Naef M, Schmid RA, et al. Prognostic role of FGFR1 amplification in early-stage non-small cell lung cancer. Br J Cancer. 2014;110:2914–22. 62. Armstrong K, Ahmad I, Kalna G, Tan SS, Edwards J, Robson CN, et al. Upregulated FGFR1 expression is associated with the transition of hormone- naive to castrate-resistant prostate cancer. Br J Cancer. 2011;105:1362–9. 63. Tomlinson DC, Lamont FR, Shnyder SD, Knowles MA. Fibroblast growth factor receptor 1 promotes proliferation and survival via activation of the mitogen-activated protein kinase pathway in bladder cancer. Cancer Res. 2009;69:4613–20. 64. Gorringe KL, Jacobs S, Thompson ER, Sridhar A, Qiu W, Choong DYH, et al. High-resolution single nucleotide polymorphism array analysis of epithelial ovarian cancer reveals numerous microdeletions and amplifications. Clin Cancer Res. 2007;13:4731–9. 65. Göke F, Bode M, Franzen A, Kirsten R, Goltz D, Göke A, et al. Fibroblast growth factor receptor 1 amplification is a common event in squamous cell carcinoma of the head and neck. Mod Pathol. 2013;26:1298–306. 66. Clauditz TS, Böttcher A, Hanken H, Borgmann K, Sauter G, Wilczak W, et al. Prevalence of fibroblast growth factor receptor 1 (FGFR1) amplification in squamous cell carcinomas of the head and neck. J Cancer Res Clin Oncol. 2018;144:53–61. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Genome Medicine Springer Journals

Loading next page...
 
/lp/springer_journal/integrated-time-course-omics-analysis-distinguishes-immediate-wm6NI7o81Y
Publisher
BioMed Central
Copyright
Copyright © 2018 by The Author(s).
Subject
Biomedicine; Human Genetics; Metabolomics; Bioinformatics; Medicine/Public Health, general; Cancer Research; Systems Biology
eISSN
1756-994X
D.O.I.
10.1186/s13073-018-0545-2
Publisher site
See Article on Publisher Site

Abstract

Background: Targeted therapies specifically act by blocking the activity of proteins that are encoded by genes critical for tumorigenesis. However, most cancers acquire resistance and long-term disease remission is rarely observed. Understanding the time course of molecular changes responsible for the development of acquired resistance could enable optimization of patients’ treatment options. Clinically, acquired therapeutic resistance can only be studied at a single time point in resistant tumors. Methods: To determine the dynamics of these molecular changes, we obtained high throughput omics data (RNA-sequencing and DNA methylation) weekly during the development of cetuximab resistance in a head and neck cancer in vitro model. The CoGAPS unsupervised algorithm was used to determine the dynamics of the molecular changes associated with resistance during the time course of resistance development. Results: CoGAPS was used to quantify the evolving transcriptional and epigenetic changes. Applying a PatternMarker statistic to the results from CoGAPS enabled novel heatmap-based visualization of the dynamics in these time course omics data. We demonstrate that transcriptional changes result from immediate therapeutic response or resistance, whereas epigenetic alterations only occur with resistance. Integrated analysis demonstrates delayed onset of changes in DNA methylation relative to transcription, suggesting that resistance is stabilized epigenetically. Conclusions: Genes with epigenetic alterations associated with resistance that have concordant expression changes are hypothesized to stabilize the resistant phenotype. These genes include FGFR1, which was associated with EGFR inhibitors resistance previously. Thus, integrated omics analysis distinguishes the timing of molecular drivers of resistance. This understanding of the time course progression of molecular changes in acquired resistance is important for the development of alternative treatment strategies that would introduce appropriate selection of new drugs to treat cancer before the resistant phenotype develops. Keywords: Acquired resistance, Precision medicine, Data integration, Time course analysis, Genomics, Epigenetics * Correspondence: Christine.Chung@moffitt.org; ejfertig@jhmi.edu Genevieve Stein-O’Brien, Luciane T. Kagohara and Sijia Li contributed equally to this work. Department of Oncology, Sidney Kimmel Comprehensive Cancer Center, Johns Hopkins University, Baltimore, MD, USA Full list of author information is available at the end of the article © The Author(s). 2018 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated. Stein-O’Brien et al. Genome Medicine (2018) 10:37 Page 2 of 22 Background untreated controls, enabling the distinction of the molecu- Recent advances to identification of gene regulation in lar mechanisms associated with acquired resistance from cancer have enabled the selection of targeted therapies those that would occur due to the long-term culturing to inhibit specific regulators of oncogenic signaling path- over the time period that resistance develops. ways essential for tumor development and maintenance Along with novel time course datasets, inferring the [1]. These therapies prolong survival but are not cura- specific and targetable signaling changes that drive tive, since most patients develop acquired resistance therapeutic resistance also requires new bioinformatics within the first few years of treatment [2]. Although a pipelines to analyze and visualize these data. The bio- wide variety of molecular alterations that confer resist- informatics pipelines must integrate genetic, epigenetic, ance to the treatment have been described, the mecha- and transcriptional changes from multiple-high through- nisms and timing of their evolution are still poorly put platforms to infer the complex gene regulatory characterized [3, 4]. As serial biopsies along the treat- mechanisms that are responsible for acquired resistance. ment period are impractical due to the invasiveness and Current supervised bioinformatics algorithms that find high costs of the procedure, the molecular alterations as- time course patterns in genomic data adjust linear sociated with acquired resistance are only known when models to correlate molecular profiles with known tem- resistance has already developed and little is known poral patterns [12–15]. Many unknown variables such as about what changes occur at earlier or later time points culture conditions, immediate response to cetuximab, during the targeted therapy. The lack of adequate in and adaptive changes may have confounding effects on vitro and in vivo time course datasets makes it challen- known covariates of therapeutic response such as growth ging to delineate the two predominant hypotheses for rates, colony size, or apoptosis rates. Unsupervised bio- how therapeutic resistance develops: (1) the presence of informatics algorithms learn the dynamics directly from small populations of resistant cells that will survive the the high-throughput data, and therefore do not require a treatment and repopulate the tumor; or (2) the develop- priori knowledge of the complex dynamics associated ment of de novo resistance mechanisms by the tumor with therapeutic response. Some unsupervised algo- cells [4, 5]. Characterization of the dynamics of genomic rithms [16–21] seek breaking points of coherent, regula- alterations induced during acquired resistance can iden- tory relationships to infer the dynamics along pathways. tify targetable oncogenic drivers and determine the best Many of these algorithms trace individual phenotypes or time point to introduce alternative therapeutic strategies individual genomics platforms. Their ability to determine to avoid resistance establishment [6]. drivers of gene expression associated with acquired re- Epidermal growth factor receptor (EGFR) inhibitors sistance from time course data in multiple experimental represent a common class of targeted therapeutics. conditions and multiple genomics data modalities is Cetuximab, a monoclonal antibody against EGFR, is emerging [22]. Further extensions are needed to contrast FDA approved for the treatment of metastatic colorectal the dynamics of signaling response to therapy to the dy- cancer and head and neck squamous cell carcinoma namics of control conditions to distinguish the specific (HNSCC) [7]. As with other targeted therapies, stable molecular processes that are unique to resistance. response is not observed for a long period and virtually Matrix factorization algorithms are unsupervised and all patients invariably develop acquired resistance [8]. can distinguish the relative molecular changes in each Recent advances in the establishment of in vitro models experimental condition over time without requiring of acquired cetuximab resistance [9] provide a unique prior knowledge of gene regulation. We have found that opportunity to study the time course of genetic events Bayesian, non-negative matrix factorization algorithms resulting in acquired resistance. Cell lines chronically ex- such as Coordinate Gene Activity in Pattern Sets posed to the targeted agent develop resistance and can (CoGAPS) [23] can extend beyond clustering to robustly be sequentially collected during the course of treatment quantify the dynamics and infer the gene regulatory net- to evaluate the progressive molecular changes. Previous works directly from the input time course data [24]. The studies to assess the mechanisms of acquired cetuximab CoGAPS error model can also be modified to enable data- resistance have been limited to comparing the genomic driven inference in distinct molecular platforms for infer- profile of the parental intrinsic sensitive cell line to ence of epigenetic regulation of gene expression [25]. stable clones with acquired resistance [9–11]. Therefore, In this study, we developed a new bioinformatics ana- these studies fail to capture the dynamics of acquired mo- lysis pipeline for integrated analysis of gene expression lecular alterations during the evolution of therapeutic re- and DNA methylation changes that occur during the sistance. The development of in vitro time course data to time course progression of resistance to targeted therap- determine the molecular drivers of therapeutic resistance ies using CoGAPS. Genes uniquely associated with these is crucial. These experimental systems have the further ad- changes were selected using a PatternMarker statistic vantage that time course data can also be generated for [26] to enable novel visualization of molecular Stein-O’Brien et al. Genome Medicine (2018) 10:37 Page 3 of 22 alterations dynamics inferred with CoGAPS. In order to three days for 11 weeks (generations G1 to G11). On the benchmark our new bioinformatics pipeline, we used an eighth day, cells were harvested. Sixty thousand cells in vitro HNSCC cell line model to induce resistance and were replated for another week of treatment with cetuxi- measure the molecular changes using high-throughput mab and the remaining cells were separately collected assays while the resistant phenotype developed. Gene ex- for: (1) RNA isolation (gene expression analysis); (2) pression and DNA methylation changes were screened DNA isolation (DNA methylation analysis); (3) prolifera- weekly while acquired cetuximab resistance was induced tion assay; and (4) storage for future use. All steps were in SCC25 cell line (intrinsic sensitive to cetuximab) and repeated for a total of 11 weeks. In parallel with the compared to the status of the untreated controls at the cetuximab treated cells, we generated controls that re- same culturing time point. CoGAPS [26] inferred spe- ceived the same correspondent volume of phosphate cific patterns of expression and DNA methylation that buffered saline (PBS). Cells were plated in several repli- are associated with the gradual establishment of ac- cates each time at the same initial density. The replicates quired cetuximab resistance. The onset of methylation were then harvested and pooled to provide enough cells changes associated with resistance is temporally delayed for genetic, epigenetic, and proliferation assays. To relative to expression changes suggesting that epigenetic achieve adequate final cell confluence and number of alterations stabilize the transcriptional changes relevant cells for the experimental analysis of each generation, to the resistant phenotype. This analysis found anti- cetuximab- and PBS-treated cells were plated in differ- correlated changes between DNA methylation and gene ent flask sizes. Cells treated with cetuximab were plated expression in FGFR1 during acquired therapeutic resist- in multiple T75 (75cm ) flasks (60,000 cells/flask) that ance. Upregulation of FGFR1 has previously been associ- were combined on the eighth day. PBS-treated cells were ated as a mechanism of acquired cetuximab resistance in plated in a single T175 (175cm ) flask (60,000 cells). HNSCC [27–29]. The identification of a canonical marker This design was selected considering the growth of resistance to EGFR inhibitors in this present study cor- inhibition of the earliest cetuximab generations and to roborates the efficacy of our experimental model and control confluence of the PBS controls at the collection analytical algorithm to identify mechanisms of resistance. time (Additional file 1: Figure S1). To our knowledge, this is the first demonstration of the anti-correlation between FGFR1 methylation and expres- Cell proliferation and colony formation assays sion suggestive of its epigenetic regulation in acquired re- Cell proliferation events were measured using the Click-iT sistance to cetuximab. Thus, this pipeline can identify Plus EdU Flow Cytometry Assay Kit Alexa Fluor 488 Pico- mechanisms of gene regulation in acquired resistance lyl Azide (Life Technologies, Carlsbad, CA, USA) accord- from high-throughput, multi-platform time course data. ing to manufacturer’s instructions. The cetuximab The resulting bioinformatics pipeline is poised to infer the generations were considered resistant when the frequency dynamics of acquired resistance from emerging time of proliferating cells was higher than in the PBS control course data with other cancer types and therapeutics. generations. Proliferation curves were generated using lo- cally weighted polynomial regression (lowess) in R. Methods Anchorage-independent growth assay was used to fur- Cell lines and materials ther confirm the development of resistance. The parental SCC25 cells were purchased from American Type SCC25 and the late G10 resistant cells were treated with Culture Collection (ATCC). Cells were cultured in different concentrations of cetuximab 10 nM, 100 nM, Dulbecco’s Modified Eagle’s medium and Ham’s F12 and 1000 nM. Number of colonies was compared to the medium supplemented with 400 ng/mL hydrocortisone same cells treated with PBS. Colony formation assay in and 10% fetal bovine serum and incubated at 37 °C and Matrigel (BD Biosciences, Franklin Lakes, NJ, USA) was 5% carbon dioxide. The parental cell line SCC25 and the performed as described previously [30]. late cetuximab and PBS generation 10 were authenti- cated using short tandem repeat (STR) analysis kit Stable SCC25 cetuximab resistant single clones (CTXR clones) PowerPlex16HS (Promega, Madison, WI, USA) through Resistance to cetuximab was induced in an independent the Johns Hopkins University Genetic Resources Core passage of SCC25 cells. After resistance was confirmed, Facility. Cetuximab (Lilly, Indianapolis, IN, USA) was single cells were isolated and grown separately to generate purchased from the Johns Hopkins Pharmacy. the isogenic resistant single cell clones (CTXR). In total, 11 CTXR clones were maintained in culture without Induction of cetuximab resistance and time course sample addition of cetuximab. With the exception of one clone collection (CTXR6), all CTXR clones presented substantial survival The HNSCC cell line SCC25 (intrinsically sensitive to advantage compared to the parental SCC25, as reported cetuximab) was treated with 100 nM cetuximab every by Cheng et al. [31]. Each of these clones was Stein-O’Brien et al. Genome Medicine (2018) 10:37 Page 4 of 22 authenticated using STR analysis kit GenePrint 10 HumanMethylation450 BeadChip platform (Illumina) at (Promega) through the JHU-GRCF, as previously the JHMI Sidney Kimmel Cancer Center Microarray published [31]. Core Facility. Briefly, DNA quality was assessed using Proliferation assay was performed to confirm cetuximab the PicoGreen DNA Kit (Life Technologies) and 400 ng resistance in the CTXR clones compared to the parental of genomic DNA was bisulfite converted using the EZ SCC25. A total of 1000 cells were seeded in 96-well plates DNA Methylation Kit (Zymo Research, Irvine, CA, in quadruplicate for each condition. PBS or cetuximab USA) following manufacturer’s recommendations. A (10 nM, 100 nM or 1000 nM) was added after 24 and total volume of 4 μL of bisulfite-converted DNA was 72 h and cells were maintained in culture for seven days. denatured, neutralized, amplified, and fragmented ac- AlamarBlue reagent (Invitrogen, Carlsbad, CA, USA) at a cording to the manufacturer’s instructions. Finally, 12 μL 10% final concentration was incubated for 2 h and fluores- of each sample were hybridized to the array chip cence was measured according to the manufacturer’s rec- followed by primer-extension and staining steps. Chips ommendations (545 nm excitation, 590 nm emission). were image-processed in the Illumina iScan system. Data Resistance in the CTXR clones was confirmed when the from the resulting iDat files were normalized with fun- proliferation rates were higher than in the PBS-treated norm implemented in the R/Bioconductor package minfi SCC25 cells. (version 1.16.1) [35]. Methylation status of each CpG site was computed from the signal intensity in the methyl- RNA-sequencing (RNA-seq) and data normalization ated probe (M) and unmethylated probe (U)asa β value RNA isolation and sequencing were performed for the as follows: parental SCC25 cells (G0) and each of the cetuximab and PBS generations (G1 to G11) and the CTXR clones M β ¼ : at the Johns Hopkins Medical Institutions (JHMI) Deep M þ U Sequencing & Microarray Core Facility. RNA-seq was also performed for two additional technical replicates of Annotations of the 450K probes to the human genome parental SCC25 cell line to distinguish technical variabil- (hg19) were obtained from the R/Bioconductor package ity in the cell line from acquired resistance mechanisms. FDb.InfiniumMethylation.hg19 (version 2.2.0). Probes on Total RNA was isolated from 1 × 10 cells using the sex chromosomes or annotated to single nucleotide AllPrep DNA/RNA Mini Kit (Qiagen, Hilden, Germany) polymorphisms were filtered from analysis. The CpG is- following manufacturer’s instructions. The RNA land probe located closest to the transcription start site concentration was determined by the spectrophotometer was selected for each gene. Genes with CpG island Nanodrop (Thermo Fisher Scientific, Waltham, MA, probes < 200 bp from the transcription start site were USA) and quality was assessed using the 2100 retained to limit analysis to CpG island promoter probes Bioanalyzer (Agilent, Santa Clara, CA, USA) system. An for each gene. Probes were said to be unmethylated for RNA Integrity Number (RIN) of 7.0 was considered as β < 0.1 and methylated for β > 0.3 based upon thresholds the minimum to be used in the subsequent steps for defined in TCGA analyses [34]. All DNA methylation RNA-seq. Library preparation was performed using the data from this study are available from GEO (GSE98813) TrueSeq Stranded Total RNAseq Poly A1 Gold Kit (Illu- as part of SuperSeries GSE98815. mina, San Diego, CA, USA), according to manufacturer’s recommendations, followed by messenger RNA (mRNA) Hierarchical clustering and CoGAPS analysis enrichment using poly(A) enrichment for ribosomal The following filtering criterion for genes from the profil- RNA (rRNA) removal. Sequencing was performed using ing of the time course data from generations of cetuximab the HiSeq platform (Illumina) for 2 × 100 bp sequen- treated cells was used. Genes from RNA-seq data were se- cing. Reads were aligned to hg19 with MapSplice [32] lected if they had log fold change > 1 between any two and gene expression counts were quantified with RSEM time points of the same condition and < 2 between the [33]. Gene counts were upper-quartile normalized and replicate control samples at time zero (5940 genes). CpG log transformed for analysis following the RSEM v2 island promoter probes for each gene were retained if the pipeline used to normalize TCGA RNA-seq data [34]. gene switched from unmethylated (β <0.1) to methylated All RNA-seq data for the cell line mode in this study are (β > 0.3) in any two samples of the time course (1087 available from GEO (GSE98812) as part of SuperSeries genes). We used the union of the sets of genes retained GSE98815. from these filtering criteria on either data platform for analysis, leaving a total of 6445 genes in RNA-seq and DNA methylation hybridization array and normalization 4703 in DNA methylation. Genome-wide DNA methylation analysis was performed Hierarchical clustering analysis was performed with on the same samples as RNA-seq using the Infinium Pearson correlation dissimilarities between genes and Stein-O’Brien et al. Genome Medicine (2018) 10:37 Page 5 of 22 samples on all retained genes. CoGAPS analysis was per- meta-pathway were derived for each pattern using the formed on both log transformed RNA-seq data and PatternMarkers statistics [26]. Comparisons between DNA methylation β values, independently using the R/ DNA methylation and gene expression values for Bioconductor package CoGAPS [23] (version 2.9.2). PatternMarkerGenes or from CoGAPS patterns and am- CoGAPS decomposes a matrix of data D according to plitudes were computed with Pearson correlation. the model Gene set analysis of cetuximab resistance signatures, the EGFR network, and pathways D ∼N A P ; Σ ; i; j i;k k; j i; j Gene set activity was estimated with the gene set statis- k¼1 tic implemented in calcCoGAPSStat of the CoGAPS R/ Bioconductor package [23]. Analyses were performed on where N represents a univariate normal distribution, three gene sets: resistance signatures, gene targets of matrices A and P are learned from the data for a speci- transcription factors in the EGFR network, and canon- fied number of dimensions P, Σ is an estimate of the i, j ical pathways. Resistance signatures were defined based standard deviation of each row and column of the data on previous literature. Specifically, in a previous study, matrix D, and i represents each gene and j each sample. CoGAPS learned a meta-pathway from gene expression In this decomposition, each row of the pattern matrix P Val12D data corresponding to overexpression of the HRAS quantifies the relative association of each sample with a in the HaCaT cell line model. That study associated the continuous vector of relative gene expression changes in CoGAPS HaCaT-HRAS meta-pathway with gene expres- the corresponding column of A. That is each row of P sion changes in acquired cetuximab resistance in the provides the relative magnitude across samples are called HNSCC cell line UMSCC1 [23]. In the current study, we patterns and quantify the separation of distinct experi- applied the PatternMarkers statistic [26] to the previ- mental conditions. These relative gene weights in the ously published CoGAPS analysis of these data to derive columns of A represent the degree to which each gene is a gene set from this meta-pathway called HACAT_ associated with an inferred pattern and are called meta- HRAS_CETUXIMAB_RESISTANCE or HACAT_RESIS pathways. Together, these matrices provide a low- TANCE. In addition, we searched MSigDB [37] (version dimensional representation that reconstructs the signal 5.2) for all gene sets associated with resistance to EGFR of the input genomics data. A single gene may have inhibition. In this search, we found the gene sets COL non-zero magnitude in several distinct gene sets, repre- DREN_GEFITINIB_RESISTANCE_DN and COLDREN_ senting the fact that a single gene can have distinct roles GEFITINIB_RESISTANCE_UP representing resistance in different biological processes (such as immediate to the EGFR inhibitor gefitinib in non-small-cell lung can- therapeutic response and acquired resistance). A recently cer (NSCLC) cell lines [38]. Gene sets of transcription fac- developed PatternMarker statistic [26] selects the genes tor targets were obtained from experimentally validated that are unique to each of the inferred patterns and targets annotated in the TRANSFAC [39] professional therefore represent biomarkers unique to the corre- database (version 2014.1). Canonical pathways were ob- sponding biological process. tained from the C2 set of MSigDB [37](version6.1). In the CoGAPS analysis of the data in this study, the standard deviation of the expression data was 10% of the Sources and analysis of additional in vitro and human signal with a minimum of 0.5. The standard deviation of tumor genomics data DNA methylation data under the assumption that β Genomics analyses of TCGA were performed on level 3 values follow a beta distribution is RNA-seq and DNA methylation data from the 243 HPV- vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi negative HNSCC samples from the freeze set for publica- β 1−β i; j i; j tion [36]. DNA methylation data were analyzed for the Σ ¼ : i; j M þ U þ 1 same CpG island promoter probes obtained in the cell line i; j i; j studies. Pearson correlation coefficients were computed in CoGAPS was run for a range of 2–10 dimensions P R to associate different molecular profiles. for expression and 2–5 for DNA methylation. Robust- Additional analysis was performed on Affymetrix Hu- ness analysis with ClutrFree [36] determined that the man Genome U133 plus 2.0 GeneChip arrays for the optimal number of dimensions P for expression was 5. SCC1/1CC8 isogenic cetuximab sensitive and resistant DNA methylation was run in four parallel sets using cell line pair described previously (GEO GSE21483 [30]). GWCoGAPS [26]. In DNA methylation, the maximum Additional gene expression data from SCC25 generated number of patterns that modeled resistance mechanisms from the same platform in the same lab was also used over and above technical variation in replicate samples for analysis, using fRMA for normalization [40] to con- of SCC25 was three. Gene sets representative of the trol for batch effects as described previously [41]. Stein-O’Brien et al. Genome Medicine (2018) 10:37 Page 6 of 22 Analysis was also performed on gene expression data 2000). It was prospectively planned to perform transla- measured with Illumina HumanHT-12 WG-DASL V4.0 tional research and patients gave their informed consent R2 expression beadchip arrays on pre-treatment samples for repeated biopsies. from patients subsequently treated with cetuximab from Bossi et al. [42], using expression normalization and Results progression-free survival groups as described in the Prolonged exposure to cetuximab induces resistance study. Data were obtained from the GEO GSE65021 Cetuximab resistance was induced by treating the series matrix file. SCC25 cells for a period of 11 weeks (CTX-G1 to – DNA samples from eight human tumor surgical speci- G11). SCC25 cells treated with PBS were used as time-- men post cetuximab treatment from the sample cohort matched controls (PBS-G1 to –G11). Response to cetux- in Schmitz et al. [43] were obtained for methylation pro- imab was determined by comparing the proliferation filing. Specifically, for each tumor one FFPE slide was rates between CTX and PBS generations. Proliferation of stained with hematoxylin and eosin and tumor burden the PBS generations is stable throughout the 11 weeks was evaluated. When the tumor content was < 50%, the (G1 to G11, Fig. 1a). Conversely, proliferation of the adjacent unstained FFPE slides were macrodissected in CTX generations progressively increases over each week order to enrich the tumor burden. A double code was (Fig. 1a). Relative to the untreated controls, the growth assigned to each sample by Biorepository. DNA was then of the treated cells is initially (CTX-G1) inhibited until extracted from two unstained slides using the QIAamp CTX-G3. Starting at CTX-G4, the absolute proliferation DNA FFPE Tissue kit. Briefly, slides were dipped into a values are equal at this week, but the fit to the data xylene bath until paraffin was melted. Then slides were across all time points suggests that the cells become re- washed with ethanol 100% and tissue was harvested for sistant to the anti-proliferative effects of cetuximab and extraction with QIAGEN affinity columns. The extracted gain stable growth advantages compared to the un- DNA was quantified by NanoDrop spectrophotometer. treated controls (CTX-G8 to –G11). DNA methylation was measured with the Illumina Comparison of proliferation rates between generations MethylationEPIC BeadChip (850K) array. Array data of CTX-treated cells relative to generations of cells were normalized with the NOOB method [44] and con- treated with PBS enabled us to conclude that cell growth verted to virtual 450K arrays using the R/Bioconductor advantages arise from chronic cetuximab treatment and package minifi version 1.22.1 and are available from are associated with resistance rather than prolonged cell GEO SuperSeries (GSE110995). Two samples had DNA culturing. We mirrored the changes in proliferation rates content < 250 ng and clustered separately from the with clinical responses seen in HNSCC tumors treated remaining samples. These samples were filtered as low with cetuximab (Fig. 1b). The lower growth rates in quality and excluded from the analysis, leaving six total CTX-G1 to -G3 may be an equivalent to the initial ef- tumor samples with DNA methylation data. Probes se- fects of the clinical treatment when the cancer cells are lected for the in vitro Illumina 450K DNA methylation sensitive to cetuximab and reduction of the tumor size is data were used for subsequent analyses. Gene expression observable. Even with gain in cell proliferation at CTX- data from biopsy samples before cetuximab treatment G3 and -G4, our model still corresponds to response to and surgical samples after cetuximab treatment were ob- the treatment since the treated generations are not tained from the previous Schmitz et al. [43] study and growing more than the controls (clinical stable tumor normalized as described previously [41] and available size). Finally, from CTX-G4 the higher proliferation even from GEO SuperSeries (GSE110996). with cetuximab treatment is a representation of acquired We performed t-tests and projections in R on the resistance noticeable in the HNSCC patients as tumor probe that had the highest standard deviation of expres- recurrence or increase in tumor size. sion values for each gene. CoGAPS signatures were also Higher proliferation in treated than in untreated cells projected into these gene expression data using the starts at CTX-G4 and we established this time point to methods described in Fertig et al. [23] with the ProjectR call as the moment at which cetuximab resistance is sta- package version 0.99.15 available from Github (https:// bly acquired and all subsequent time points continue to github.com/genesofeve/projectR). We also performed t- develop acquired stable cetuximab resistance. To con- tests in R to compare the long- (LPFS) and short-term firm this hypothesis, we evaluated the ability of the re- progression-free survival (SPFS) groups based on the sistant CTX-G10 to anchorage-independent growth. values obtained from this projection. Even under different concentrations of cetuximab, HNSCC samples and patient information collection CTX-G10 presents enhanced anchorage-independent were approved by the Independent Ethics Committee growth compared to the parental SCC25 (G0) (two-way and the Belgian Health Authorities and conducted in ac- ANOVA with multiple comparisons p value < 0.01 for cordance with the Declaration of Helsinki (October each concentration, Additional file 1: Figure S2), Stein-O’Brien et al. Genome Medicine (2018) 10:37 Page 7 of 22 a Proliferation assay Generations (weeks of treatment) b Clinical evolution Tumor detection Treatment with cetuximab Tumor reduction Stable tumor size Tumor growth / Disease recurrence RESPONSE TO CETUXIMAB ACQUIRED RESISTANCE Unaccessible molecular changes Fig. 1 In vitro time course reflects clinical evolution of cetuximab response and evolution of acquired resistance. Intrinsic cetuximab-sensitive HNSCC cell line SCC25 was treated with cetuximab (red) or PBS (black) for 11 generations to develop acquired resistance. Proliferation assay (flow) of cetuximab treatment (red line) and PBS-treated cells (black line) measured cetuximab response for all SCC25 generations. While proliferation of the PBS generations was stable throughout the 11 weeks, proliferation of the CTX generations progressively increased over each week. Relative to the untreated controls, the growth of the treated cells was initially (CTX-G1) inhibited until CTX-G3. Starting at CTX-G4, the cells became resistant to the anti-proliferative effects of cetuximab and gained stable growth advantages compared to the untreated controls demonstrating the stabilization of cetuximab resistance proliferation rates than the PBS controls (shown in in later generations. Fig. 1a). The two groups of samples resistant to cetuximab are represented by a progressive increase Treatment vs control gene expression changes governs in proliferation that is more significant than in the clustering and immediate therapeutic response is untreated controls (weeks 4 to 8) and by the stabilization confounded with changes from acquired resistance in the proliferation rates (weeks 9 to 11), but still higher RNA-seq data for the parental SCC25 cell line (G0) and than in thePBS generations. Theexpressionchanges at from each generation of CTX- and PBS-treated cells the distinct time points during development of acquired were collected to characterize the gene expression resistance are shared among numerous genes. Although changes occurring as cells acquired cetuximab resist- the clustering was able to separate cetuximab- from ance. Gene expression changes between treated (cetuxi- PBS-treated cells, it was not possible to discriminate mab) and untreated (PBS) cells and over generations of the alterations related to an immediate therapeutic re- treated cells are apparent in time-ordered RNA-seq data sponse (not relevant to the resistant phenotype) from (Fig. 2a). Additional clustering analysis of the samples resistance-specific gene expression changes. accounting for the treatment time point (generations/ Similar separation of the sensitive and resistant phases columns) (Additional file 1: Figure S3) distinguish three of cetuximab response is observed in clustering analysis clusters of samples: those with cetuximab sensitivity of gene signatures previously described in HNSCC and (CTX-G1 to CTX-G3); those with early cetuximab re- NSCLC cell line models resistant to cetuximab or gefi- sistance (CTX-G4 to CTX-G8); and those with late or tinib (anti-EGFR small molecule), respectively [38, 41] stable cetuximab resistance (CTX-G9 to CTX-G11). The (Additional file 1: Figure S4). For these genes, changes group of cetuximab sensitive samples corresponds to the during early phases of resistance clusters for CTX-G4 to time points at which the CTX generations present lower CTX-G6 as distinct from later generations CTX-G7 to Cell proliferation (%) 10 15 20 25 30 35 40 Stein-O’Brien et al. Genome Medicine (2018) 10:37 Page 8 of 22 Clustering of gene expression data CoGAPS analysis of gene expression data PBS Cetuximab Expression of PatternMarker genes bc CoGAPS expression patterns for CoGAPS expression patterns (sample weights v. time) PBS Cetuximab 0 24 68 10 Color Key Row Z-score (log2 RSEM) Weeks of treatment PBS −4 0 4 Cetuximab weeks of treatment Fig. 2 CoGAPS analysis identifies signatures of resistance to EGFR inhibitors and separate resistant and control generations. a Heatmap of gene expression values in 11 generations of SCC25 cells treated with 100 nM of cetuximab (red columns) to acquire resistance and with PBS as control (black columns). b Heatmap of gene expression values for PatternMarker genes identified with CoGAPS analysis of gene expression data from 11 generations of SCC25 cells treated with PBS as control (black columns) and with 100 nM of cetuximab (red columns) to acquire resistance. Rows are colored according to which CoGAPS pattern the PatternMarker statistic assigned each gene and sorted by the PatternMarker statistic. c CoGAPS patterns inferred from gene expression data over generations of PBS control (black lines) or treatment with 100 nM of cetuximab (red lines) CTX-G11. Nevertheless, these signatures also cluster gene can have high amplitude in multiple patterns, mod- samples with gene expression changes at early phases eling multiple regulation of genes but complicating (CTX-G1 to CTX-G3) as distinct from samples from visualization of the inferred patterns. A recently devel- PBS-treated generations. However, these analyses were oped PatternMarker statistic defines genes that are insufficient to quantify the relative dynamics of genes as- uniquely associated with each of these patterns. Limiting sociated with immediate response to therapy or subse- the heatmap to these genes enables visualization of the quent acquired resistance. dynamics of gene expression changes in our time course dataset (Additional file 1: Figure S5). CoGAPS analysis of gene expression distinguishes patterns In this heatmap of CoGAPS PatternMarker genes, we of acquired resistance from immediate therapeutic response observe that only three patterns (EP1, EP2, and EP3) dis- To define gene expression signatures for treatment effect tinguish the experimental conditions (cetuximab vs PBS) and cetuximab resistance, we applied the CoGAPS [26] (top three patterns on Additional file 1: Figure S5). The Bayesian matrix factorization algorithm to the time other two patterns, EP4 and EP5, represent changes in course gene expression data. CoGAPS decomposes the gene expression from the parental cell lines and subse- input data into two matrices: a pattern matrix with rela- quent generations or an expression pattern that is con- tive sample weights along rows and an amplitude matrix stant and corresponds to signature of highly expressed with relative gene weights along columns. Each row of genes (lower two patterns on Additional file 1: Figure the pattern matrix quantifies the extent of transcrip- S5), respectively. We note that highly expressed genes tional changes within the genes in the corresponding associated with EP5 may also have dynamic changes due column of the amplitude matrix and provides a low di- to treatment and are filtered in the PatternMarker ana- mensional representation of the biological process in lysis of all patterns in Additional file 1: Figure S5. EP4 that data. In this analysis, we identified five CoGAPS represents expression changes between treated cells and patterns (Expression Patterns [EP]) (Additional file 1: the parental cell line, which have a technical effect on Figure S5) in the time course gene expression dataset. A gene expression. Notably, even the exclusion the flat 11 Stein-O’Brien et al. Genome Medicine (2018) 10:37 Page 9 of 22 pattern for highly expressed genes (EP5) still retains expression signatures associated consistently with the highly expressed genes in the signature. Retaining the development of cetuximab resistance. Gene set statistics technical pattern (EP4) in this calculation filters genes of transcription factor targets of EGFR on CoGAPS gene with expression changes from technical artifacts in the weights are significantly downregulated in this acquired experimental conditions. The resulting set of Pattern- resistance pattern (Additional file 1: Figure S6). One strik- Marker genes for EP1–EP3 enable visualization of the ing exception is c-Myc, which trends with acquired resist- expression changes dynamics that are associated with ance (p value of 0.06), consistent with the role of this cetuximab response (Fig. 2b) and allow the definition of transcription factor in cellular growth. Resistance signa- a gene signature associated with that response ture COLDREN_GEFITINIB_RESISTANCE_DN is sig- (Additional file 2: Table S1). nificantly downregulated in EP2 (p value of 0.04). There Similar to the separation seen with clustering are 32 statistically significant canonical pathways associ- (Additional file 1: Figure S5), the CoGAPS pattern EP1 ated with this pattern, including notably telomerase, PI3K, distinguishes cetuximab from PBS at every generation and cell cycle pathways (Additional file 3: Table S2). (Fig. 2b and c, top). These genes present an immediate The third CoGAPS expression pattern (EP3) repre- transcriptional upregulation in response to cetuximab sents a gradual repression of gene expression with treatment. Gene set analysis to determine the function cetuximab treatment (Fig. 2b and c, bottom). This of CoGAPS patterns was performed with an enrichment expression pattern trends to significant enrichment in the analysis on all gene weights in the amplitude matrix ob- COLDREN_GEFITINIB_RESISTANCE_DN resistance tained from the CoGAPS analysis. By performing the signature (Additional file 1:Figure S6,one-sided p value 0. analysis on gene weights and not only the PatternMarker 12) and downregulated in the HACAT_HRAS_CETUXI- genes, as shown in Fig. 2b, we account for multiple regu- MAB_RESISTANCE resistance signature (Additional file 1: lation of genes in pathways. Specifically, we performed Figure S6, one-sided p value 0.09). This confirms that EP3 gene set analysis on published resistance signatures [38, is associated with repression of gene expression during ac- 41], transcription factors previously associated with the quired cetuximab resistance. There are also 29 statistically EGFR signaling network during cetuximab response in significant canonical pathways associated with this pat- HNSCC [39, 41], and canonical pathways from MSigDB tern, including cell lineage, metabolic, WNT, and GSK3 [37, 45] (Additional file 1: Figure S6; Additional file 3: pathways (Additional file 3: Table S2). Table S2). Gene set analysis confirms that published re- sistance signatures [38, 41] are significantly enriched in EP1 (Additional file 1: Figure S6; one-sided p values of 0. Changes in DNA methylation inferred with CoGAPS are 002 and 0.003 for resistance gene sets COLDREN_GEFI- associated with resistance to cetuximab, but not the TINIB_RESISTANCE_DN and HACAT_HRAS_CETUX immediate response to treatment observed in gene IMAB_RESISTANCE, respectively). However, the tran- expression scriptional changes in this pattern are not associated To determine the timing of the methylation changes as- with acquired resistance to cetuximab and even decrease sociated with acquired resistance, we also measured modestly as resistance developed. Further, enrichment DNA methylation in each cetuximab generation of by transcription factor AP-2alpha targets (TFAP2A; one- SCC25 cells and PBS controls (Fig. 3a). Application of sided p value of 0.05) confirms previous work indicating the CoGAPS matrix factorization algorithm with the that transcription by AP-2alpha is induced as an early PatternMarker statistics to the methylation data reveals feedback response to EGFR inhibition [39]. There are 84 a total of three methylation patterns (MP) (Fig. 3bc; significant canonical pathways from MSigDB, including Additional file 2: Table S1): gradual increase of DNA notably pathways associated with the immune system, methylation in controls (MP1, Fig. 3b middle); rapid extracellular matrix, ERBB4 signaling, and VEGF signal- demethylation in CTX generations starting at CTX-G4 ing (Additional file 3: Table S2). Based upon these (MP2, Fig. 3b bottom); and rapid increase in DNA findings, we concluded that EP1 is associated with im- methylation in CTX generations starting at CTX-G4 mediate response to cetuximab although it includes (MP3, Fig. 3b top). In contrast to the gene expression genes that are also associated with cetuximab resistance data, there is no immediate shift in DNA methylation in previous studies. resulting from cetuximab treatment. Gene set analysis The second CoGAPS expression pattern (EP2) quanti- was performed on canonical pathways from MSigDB fies divergence of the cetuximab treated cells from (Additional file 4: Table S3) and found 26 statistically controls at generation CTX-G4 (Fig. 2b and c, middle) significant pathways for MP1, 29 for MP2, and 27 for which is the time point that cetuximab treated cells MP3. In contrast to gene expression, the majority of ca- present significant and stable growth advantage over nonical pathways is shared by the three methylation PBS controls (Fig. 1a). Therefore, EP2 contains gene patterns and include notably the cytokine (PID-IL8- Stein-O’Brien et al. Genome Medicine (2018) 10:37 Page 10 of 22 Clustering of DNA methylation data CoGAPS analysis of DNA methylation data PBS Cetuximab DNA methylation of PatternMarker genes CoGAPS DNA methylation patterns bc for CoGAPS methylation patterns (sample weights v. time) PBS Cetuximab 1 1 weeks of treatment 0 24 68 10 Weeks of treatment Color Key Row Z-score (b) PBS Cetuximab −4 0 4 weeks of treatment Fig. 3 Dynamics of DNA methylation alterations and patterns in acquired cetuximab resistance. a Heatmap of DNA methylation values in 11 generations of SCC25 cells treated with PBS as control (black columns) and with 100 nM of cetuximab (red columns) to acquire resistance. b Heatmap of DNA methylation values for genes selected by CoGAPS DNA methylation patterns analysis in the same SCC25 cetuximab and PBS generations. c CoGAPS patterns inferred from DNA methylation data over generations of PBS control (black lines) or treatment with 100 nM of cetuximab (red lines) CXCR2 and IL8-CXCR1 pathways) and FGFR of the expression and methylation patterns previously (Reactome Signaling by FGFR3 mutants) signaling shown (Figs. 2c and 3c, respectively) enable visualization pathways. of the time point when changes between cetuximab and Comparing the CoGAPS patterns from gene expression PBS generations are significant in each pattern. These and DNA methylation reveals strong anti-correlation be- dynamics explain the discrepancy between the genes as- tween gene expression and DNA methylation in resistant sociated with each pattern and suggest that DNA methy- patterns (Additional file 1: Figure S7A). We observed that lation stabilizes the gene expression signatures crucial to the gene expression changes associated with acquired re- the maintenance of acquired cetuximab resistance. sistance occur gradually and are evident in early genera- tions (Fig. 2c). The DNA methylation is consistent in Gene expression and methylation profile of SCC25 cetuximab treatment and control PBS in DNA methyla- single-cell clones with acquired cetuximab resistance tion patterns MP2 and MP3 during early generations; fol- demonstrates cell heterogeneity lowing which, there is a rapid accumulation in DNA The little overlap between the gene expression and DNA methylation changes starting after generations CTX-G4 methylation PatternMarker genes and non-specific DNA and CTX-G5 in both MP2 and MP3 (Fig. 3c), concurrent methylation pathways may arise due to the development with the onset of the observed growth advantage over the of different resistant sub-clones with specific gene signa- PBS control (Fig. 1a). tures of acquired resistance in bulk data. In order to While the patterns themselves are anti-correlated, the address this issue and to delineate whether our pre- gene weights that define meta-pathways and are inferred sumptive drivers resulted from clonal expansion of re- in the amplitude matrix corresponding to each pattern sistant cells or from the development of new epigenetic with CoGAPS are not (Additional file 1: Figure S7B). alterations to drive resistance, we measured DNA We also observed little overlap between the PatternMar- methylation and gene expression on a panel of 11 iso- ker genes from methylation patterns and gene expres- genic stable cetuximab-resistant clones (CTXR1 to sion. Changes in DNA methylation are delayed relative CTXR11) derived from SCC25 cells in a previous study to those of gene expression in acquired cetuximab resist- [31]. Despite being derived from the parental SCC25 ance as can be noted in Fig. 4, where direct comparison cells after chronic exposure to cetuximab, the CTXR 11 Stein-O’Brien et al. Genome Medicine (2018) 10:37 Page 11 of 22 Fig. 4 DNA methylation and expression CoGAPS patterns demonstrate delayed onset of epigenetic changes in acquired resistance. CoGAPS patterns for gene expression (top) and DNA methylation (bottom) of patterns associated with acquired cetuximab resistance in SCC25 cetuximab generations (red) relative to PBS generations (black). Vertical dashed line represents time at which patterns for SCC25 generation separated from pattern for PBS generations. The timing of methylation changes distinguishing cetuximab-resistant generations was delayed in DNA methylation relative to that of gene expression clones and the time course generations display wide- analysis between DNA methylation and gene expression for spread differences. We plot both expression and methy- each of the DNA methylation PatternMarker genes (Fig. 5). lation profiles (Additional file 1: Figure S8 and Figure This analysis identified FGFR1 as one of the genes S9, respectively) among the DNA methylation Pattern- with significant anti-correlation between expression Marker genes that are anti-correlated with expression and methylation, suggesting potential epigenetic regu- and cellular morphology (Additional file 1: Figure S10) lation during cetuximab resistance acquisition. This for these clones. Significantly greater heterogeneity is finding is consistent with previous studies that associ- observed among the CTXR clones in all these genomic ate differential expression of FGFR1 with resistance to and morphological data. Figure 5 demonstrates that EGFR inhibitors, including cetuximab, in HNSCC and higher heterogeneity among single cell clones is also ob- other tumor types in vitro and in vivo [27, 46–48]. served in the epigenetically regulated PatternMarker However, none of these studies demonstrate an asso- genes from the CoGAPS analysis. These results suggest ciation between FGFR1 upregulation and demethyla- that different mechanisms of resistance may arise in the tion. Given the tight temporal regulation of the DNA same HNSCC cell line as a result of intra-heterogeneity, methylation PatternMarkers with anti-correlated expres- resulting in the detection of a wide range of expression sion and the previous work on FGFR1,we hypothesize signatures with higher or lower correlation with the that this set of genes represents epigenetic drivers of methylation profile depending on the size of each spe- acquired resistance. cific cell population. We hypothesize that epigenetically regulated genes shared along the time course patterns and resistant FGFR1 overexpression and demethylation are associated single-cell clones might implicate common mechanisms with acquired cetuximab resistance in the time course acquired during evolution of the stable resistance pheno- and in stable cetuximab-resistant clones type. To test this assumption, we also performed correl- To ascertain potential drivers of the stable cetuximab- ation analysis for each of the epigenetically regulated resistant phenotype induced by DNA methylation, we genes in our resistant set (Fig. 5) in the resistant single defined genes that are PatternMarkers [26] of the DNA cell clones and parental cell lines. Nine of the epigeneti- methylation patterns associated with stable acquired cally regulated PatternMarker genes also have significantly cetuximab resistance (MP2 and MP3). We then applied anti-correlated gene expression and DNA methylation in correlation analysis to determine genes that were epi- the stable cetuximab-resistant clones. Of these, only genetically regulated. Specifically, we performed correlation FGFR1 is demethylated and re-expressed in a CTXR Stein-O’Brien et al. Genome Medicine (2018) 10:37 Page 12 of 22 PatternMarker genes for CoGAPS methylation patterns that are anti-correlated with expression in time course data Gene expression in time course and resistant clones DNA methylation in time course and resistant clones ab PBS Cetuximab Resistant clones PBS Cetuximab Resistant clones EVC2 EVC2 PIGG PIGG ZNF721 ZNF721 RRAD RRAD FRYL FRYL ZNF141 ZNF141 GALC GALC GNPDA2 GNPDA2 MAEA MAEA MAEA MAEA ZNF229 ZNF229 LMX1B LMX1B VAMP5 VAMP5 ZNF71 ZNF71 PRIMA1 PRIMA1 IDUA IDUA ZNF569 ZNF569 KCNJ5 KCNJ5 GSTA4 GSTA4 ITGA6 ITGA6 HOPX HOPX HOPX HOPX KLF3 KLF3 SLC39A8 SLC39A8 SLC39A8 SLC39A8 HLA−DPB1 HLA−DPB1 CFLAR CFLAR ENPEP ENPEP SOCS6 SOCS6 SLC35D1 SLC35D1 LCORL LCORL SLC2A10 SLC2A10 NPNT NPNT TLL1 TLL1 RDM1 RDM1 NEUROG2 NEUROG2 ECHDC2 ECHDC2 GAB1 GAB1 ETS2 ETS2 DOCK10 DOCK10 GATM GATM RIPK3 RIPK3 HEYL HEYL SLC6A2 SLC6A2 COL2A1 COL2A1 ASRGL1 ASRGL1 MTTP MTTP MIR155HG MIR155HG TFAP2E TFAP2E ZNF483 ZNF483 ALS2CL ALS2CL GNAL GNAL FBXO15 FBXO15 DIO3 DIO3 BACH1 BACH1 SYT5 SYT5 LRP10 LRP10 ELOVL4 ELOVL4 COL24A1 COL24A1 THSD1 THSD1 KCNE3 KCNE3 FLT4 FLT4 ELMOD3 ELMOD3 TSPAN10 TSPAN10 SDCBP2 SDCBP2 TTC21A TTC21A RYR3 RYR3 ISLR2 ISLR2 CNPY1 CNPY1 WNT4 WNT4 EFCAB10 EFCAB10 IGFBPL1 IGFBPL1 IL17RB IL17RB GDF15 GDF15 GAS7 GAS7 LTBP3 LTBP3 RAD50 RAD50 KCNJ2 KCNJ2 PPARG PPARG SOX15 SOX15 FGFR1 FGFR1 HS6ST1 HS6ST1 KIFC3 KIFC3 CASP3 CASP3 IQGAP2 IQGAP2 ODC1 ODC1 ODC1 ODC1 C2orf70 C2orf70 ICAM5 ICAM5 PFKP PFKP PFKP PFKP PYCARD PYCARD CIT CIT STC2 STC2 TOX TOX TCFL5 TCFL5 FBN2 FBN2 FBN2 FBN2 FOXK2 FOXK2 SPSB4 SPSB4 GAPDHS GAPDHS LBX1 LBX1 TMEM145 TMEM145 COL8A1 COL8A1 GRB10 GRB10 CAD CAD ABCF2 ABCF2 WBP2NL WBP2NL TMEM45B TMEM45B BCR BCR TIMM50 TIMM50 Row Z-score (β) Row Z-score Weeks of treatment Weeks of treatment −4 0 4 −4 0 4 Fig. 5 Clonal heterogeneity does not reflect signature of epigenetically regulated genes observed in bulk time course analysis. a Heatmap of gene expression values for DNA methylation PatternMarker genes for acquired resistance that were anti-correlated between expression and DNA methylation (Fig. 4). Data includes 11 generations of SCC25 cells treated with PBS as control (black columns labeled PBS) and with 100 nM of cetuximab (red columns labeled cetuximab) to acquire resistance and gene expression data from independent, stable cetuximab-resistant clones in absence of cetuximab treatment (CTX-resistant clones). Gene expression heatmap on a red-yellow scale indicated in the color key. b Heatmap of DNA methylation data in conditions described in (a), on a blue-yellow scale indicated in the color key clone relative to the parental SCC25 cell line (Fig. 6a, pre and post cetuximab treatment, we further investigate remaining eight genes in Additional file 1:FigureS11). the pattern of expression and methylation of FGFR1 and In this analysis, overexpression and de-methylation of EGFR in publicly available datasets. Using gene expres- FGFR1 expression occurs in only one of the resistant sion and DNA methylation data from The Cancer clones (CTXR10). QRT-PCR of FGFR1 gene expression Genome Atlas (TCGA) for 243 HPV-negative HNSCC in CTXR10 relative to the parental cell line demon- pre-treatment samples [36], we verified that the upregu- strated a > 30-fold change (Fig. 6b). Furthermore, in the lation of EGFR and FGFR1 is not concomitant (Pearson resistant cell clone, increased levels of FGFR1 were as- correlation coefficient = − 0.06, p value = 0.33, Fig. 7a). sociated with increased levels of phospho-FGFR1 and Additionally, the negative correlation of FGFR1 gene decrease in EGFR and phospho-EGFR as assessed by expression and DNA methylation status is statistically Western blot (Fig. 6c). This clone is one of the fastest significant (Pearson correlation r of − 0.32, p value growing under cetuximab treatment (Additional file 1: < 0.0001, Fig. 7b),suggesting that FGFR1 transcription is Figure S12). This observation suggests that the bulk associated with demethylation in some HPV-negative data from the time course captured clonal outgrowth of HNSCC tumors. Since there is no treatment information a cetuximab-resistant clone with similar molecular fea- available for the TCGA dataset, we could not make as- tures (FGFR1 demethylation) to CTXR10 and that sumptions related to cetuximab resistance and whether clonal outgrowth is the dominant mechanism of resist- FGFR1 methylation is a consequence of the treatment. To ance in our model. assess this question, we collected new DNA methylation data for six HNSCC tumors after cetuximab treatment Observed FGFR1 dynamics in vitro recapitulates relationships from a cohort of HNSCC tumor samples described previ- from in vivo tumor genomicsand acquired cetuximab ously [43]. All six samples have low DNA methylation resistance values for FGFR1 (β-values in the range of 0.04–0.08, with In order to confirm that the mechanisms we found with a mean of 0.05), suggesting that the gene is unmethylated our in vitro approach are present in HNSCC samples in these samples. While there was insufficient DNA to SCC25 CTXR10 CTXR3 CTXR1 CTXR12 CTXR7 CTXR8 CTXR2 CTXR4 CTXR9 CTXR11 CTXR5 SCC25 CTXR1 CTXR12 CTXR4 CTXR11 CTXR9 CTXR5 CTXR7 CTXR8 CTXR3 CTXR10 CTXR2 Stein-O’Brien et al. Genome Medicine (2018) 10:37 Page 13 of 22 a FGFR1 gene expression vs DNA methylation in resistant clones CTXR10 CTXR1 CTXR12 CTXR11 CTXR3 CTXR4 CTXR2 CTXR8 CTXR9 CTXR7 7 CTXR5 SCC25 0.1 0.2 0.3 0.4 FGFR1 DNA methylation (beta cg20148210) b FGFR1 gene expression c Protein expression SCC25 CTXR10 40 FGFR1 pFGFR (Tyr653/654) 20 EGFR pEGFR (Tyr1068) 0 GAPDH SCC25 CTXR10 Fig. 6 Overexpression and de-methylation of FGFR1 in acquired cetuximab resistance is confirmed in stable SCC25-resistant clones. a Expression of FGFR1 relative to DNA methylation in stable cetuximab-resistant clones. b QRT-PCR of FGFR1 gene expression in CTXR10 relative to the parental cell line (> 30-fold change). c Western blot comparing FGFR1, phospho-FGFR1, EGFR, and phospho-EGFR in CTXR10 relative to the parental SCC25 cell line. In the resistant cell clone, increased levels of FGFR1 were associated with increased levels of phospho-FGFR1 and decrease in EGFR and phospho-EGFR quantify DNA methylation before treatment in these pa- Nonetheless, these findings suggest that similar molecular tients, FGFR1 gene expression increases after treatment in mechanisms may contribute to both mechanisms (intrin- four of the six tumor samples (Additional file 1:Figure sic and acquired) of cetuximab resistance in HNSCC. S13). While this cohort is small, these data and TCGA suggest that FGFR1 methylation is potentially associated CoGAPS signatures of resistance and therapeutic response with its re-expression in HNSCC tumor samples in re- replicate in an independent in vitro system and sponse to cetuximab treatment. significantly stratified patient samples with long- To determine whether FGFR1 is associated with cetux- vs short-progression-free survival imab response, we used gene expression data from To further illustrate that the results are reflective of HNSCC patients before cetuximab treatment available HNSCC in a general fashion, we evaluated the behavior from Bossi et al. [42]. After follow-up, the patients were of the two additional cell lines and human tumors in the separated in SPFS (median three months survival) or CoGAPS signatures using gene expression data available LPFS (median 19 months survival) according to time of from previously published studies. The HNSCC cell lines recurrence or metastasis development. Using this data- SCC1 and 1CC8 were chosen as the cetuximab-resistant set, we confirmed that EGFR expression in SPFS is 1CC8 was generated from the cetuximab sensitive SCC1 significantly lower than in the LPSF group (Fig. 7c) (log in a similar protocol used to establish the single cell fold change − 1.0, t-test p value 0.0003). These data clones [30]. Data from SCC25 were also included as a suggest that EGFR overexpression is associated with bet- reference. It is important to note that the treatment time ter response to cetuximab, consistent with the mechan- for the SCC1 and 1CC8 pair is on the order of hours vs ism of action of the therapeutic. The opposite was weeks, as used to generate the time course data. By pro- observed for FGFR1, with overexpression in SPFS vs jecting these data into the CoGAPS signatures, the rela- LPSF (Fig. 7d, log fold change 0.9, t-test p value 0.003). tionship between the sensitive SCC1 and resistant 1CC8 However, the Bossi et al. study [42] lacks DNA methyla- recapitulates the relationship between PBS and CTX tion data to assess whether FGFR1 was epigenetically reg- time course generations, respectively, in treatment ulated in these samples. Most patients with SPFS in this driven signatures (Fig. 8a, b, e, f). Conversely, CoGAPS dataset also had intrinsic resistance to cetuximab, instead signatures related to culture specific conditions failed to of acquired resistance studied in our in vitro model. produce meaningful differences between the lines (Fig. FGFR1 expression (log2 RSEM) Relative expression ( Ct) Stein-O’Brien et al. Genome Medicine (2018) 10:37 Page 14 of 22 FGFR1 and EGFR gene expression in HNSCC FGFR1 gene expression and DNA methylation ab and normal samples (TCGA) in HNSCC and normal samples (TCGA) Pearson r=-0.06 Pearson r=-0.32 p=0.33 p<0.0001 10 12 14 16 0.0 0.2 0.4 0.6 0.8 1.0 Normal Normal EGFR expression (log2 RSEM) FGFR1 methylation (beta - cg20148210) HNSCC HNSCC EGFR gene expression in HNSCC with FGFR1 gene expression in HNSCC with c d LPFS and SPFS LPFS and SPFS LPFS SPFS LPFS SPFS (n=14) (n=26) (n=14) (n=26) Fig. 7 FGFR1 gene expression and DNA methylation patterns were confirmed in independent HNSCC tumor sample datasets. a Scatter plot of gene expression for EGFR and FGFR1 in HPV-negative HNSCC samples from TCGA demonstrated that only a few HNSCC cases present increased levels of both genes and that there is no significant correlation between the expression of both genes concomitantly. b DNA methylation of FGFR1 was anti-correlated with expression in HPV-negative HNSCC TCGA samples, suggesting that upregulation of FGFR1 might be a result of promoter demethylation in primary tumors. c EGFR expression was significantly overexpressed in a group of HNSCC patients with LPFS relative to patients with SPFS in gene expression data from Bossi et al. d FGFR1 was significantly overexpressed in patients with SPFS relative to patients with LPFS in this same dataset. LPFS - long-progression-free survival; SPFS - short-progression-free survival 8c, d). Projections of the expression patterns from were measured for a prolonged exposure (11 weeks) to a CoGAPS into the cell line gene expression data were targeted therapeutic agent. Using our novel robust time also anti-correlated with projections of the methylation course integrated analysis approach, we characterized signatures in these same data. Gene expression data the molecular alterations during the development of ac- from HNSCC tumors from patients before their treat- quired cetuximab resistance using a HNSCC in vitro ment with cetuximab described in Bossi et al. [42] were model. Cell proliferation, gene expression, and DNA also analyzed. Projection into both the CoGAPS signa- methylation high-throughput analysis were performed tures of resistance and therapeutic response significantly weekly in equivalent cultures (cetuximab and PBS con- −5 stratified LPFS vs SPFS (p value = 5.2 × 10 and 3.1 × trol generations) as resistance developed. Over the −3 10 , respectively, (Fig. 8g, h). Conversely, projection in course of 11 weeks, it was possible to compare treated to the CoGAPS signature associated with culturing was (CTX) and untreated (PBS) cells grown under the same not significant (p value = 0.50, Fig. 8i). conditions. Applying robust bioinformatics algorithms, we discriminated changes associated with acquired re- Discussion sistance from those related to adaptive response to the Analysis of course omics data enables separation of cell culturing process and treatment. The SCC25 cell immediate treatment response and technical artifacts line model was chosen since this is one of the only two from molecular mechanisms of acquired therapeutic HNSCC cell lines previously used to generate isogenic resistance cetuximab-resistant cell lines [10]. However, this is the Although numerous short time course genomics studies first study to our knowledge to enable characterization of therapeutic response have been performed [49–51], of the transcriptional and epigenetic dynamics at the this is the first time that genetic and epigenetic changes early phases of therapeutic resistance. These phases EGFR expression (ILMN_1696521) FGFR1 expression (log2 RSEM) 678 9 10 11 12 6 7 8 9 10 11 12 13 FGFR1 expression (ILMN_2398261) FGFR1 expression (log2 RSEM) 6789 10 11 12 8 9 10 11 12 Stein-O’Brien et al. Genome Medicine (2018) 10:37 Page 15 of 22 HNSCCCellline gene expression projected into CoGAPS expression signatures 0.080 a b c 0.070 0.090 0.075 0.068 Response Response Response Res Res Res 0.085 Sen Sen Sen 0.070 Treatment Treatment Treatment 0.066 0.080 0 0 0 1 1 1 P-value P-value P-value * 2.1e-04 * 5.7e-02 * .93 0.065 0.064 0.075 # 1.5e-05 # 4.1e-03 # .83 SCC25 SCC1 1CC8 SCC25 SCC1 1CC8 SCC25 SCC1 1CC8 CellLine CellLine CellLine HNSCCCellline gene expression projected into CoGAPS metylation signatures 24.6 d e f * −4.4 24.4 −0.50 Response Response Response −4.5 Res 24.2 Res Res Sen Sen Sen −0.54 Treatment OG2 Treatment −4.6 Treatment 24.0 0 0 0 1 1 1 −0.58 P-value 23.8 P-value −4.7 P-value * .47 * 7.8e-03 * 1.9e-03 # .56 # 1.3e-03 # 1.6e-03 23.6 −4.8 SCC25 SCC1 1CC8 SCC25 SCC1 1CC8 SCC25 SCC1 1CC8 CellLine CellLine CellLine Gene expression from relapsed HNSCC patients from Bossi, et al projected into CoGAPS expression signatures g h i 0.14 0.130 0.17 0.125 p−value = 0.500447 0.13 p−value = 0.003070 0.16 p−value = 0.000052 0.120 0.15 0.12 0.115 0.14 0.11 0.110 0.13 0.10 Long Short Long Short Long Short Progression Free Survival length Progression Free Survival length Progression Free Survival length Fig. 8 CoGAPS gene signatures confirmed in independent in vitro and in vivo expression data. a–c Box plots of sample weights for HNSCC cetuximab sensitive cell lines SCC25, SCC1, and resistant 1CC8 in CoGAPS expression signatures. d–f Box plots of sample weights for HNSCC cetuximab-sensitive cell lines SCC25, SCC1, and resistant 1CC8 in CoGAPS methylation signatures. g-i Box plots of tumor gene expression profiles from relapsed HNSCC patients projected into CoGAPS expression signatures. Patient samples are stratified by length of progression-free survival Pattern.1 Pattern.1 Projected Pattern 1 Pattern.2 Pattern.2 Projected Pattern 2 Pattern.3 Pattern.3 Projected Pattern 3 Stein-O’Brien et al. Genome Medicine (2018) 10:37 Page 16 of 22 cannot be measured in patients due to the complexity of albeit to a lesser degree than lowly expressed genes. Be- early detection of resistance and obtaining repeated bi- cause the PatternMarker statistics includes genes that opsy samples. are uniquely associated with each inferred biological process, these highly expressed genes would be filtered Unsupervised time course analysis with CoGAPS quantifies, from associations with the dynamic conditions. To in- visualizes, and enables functional analysis of the dynamics clude these genes in the signatures defined in this study, of acquired therapeutic resistance across omics data EP5 was filtered from the calculation of the PatternMar- modalities ker statistic. Such filtering process is not required for Determining the dynamics of the molecular alterations heatmap-based visualization and filtering of flat patterns responsible for resistance requires integrated, time is recommend when defining gene signatures containing course bioinformatics analysis to quantify these alter- sets of genes that are most strongly associated with dy- ations. Based upon previous performance of Bayesian, namic changes. Patterns that reflect technical artifacts in non-negative matrix factorization algorithms in inferring the data, such as EP4, should be retained in the Pattern- dynamic regulatory networks for targeted therapeutics Marker analysis to limit the signatures associated with [49, 50], we selected CoGAPS [23] for analysis of gene inferred processes to retain only biologically relevant expression data from our time course experiment. genes. We note that these PatternMarker statistic are CoGAPS have already proven highly effective in relating similar to the D-scores proposed in Zhu et al. [52] and gene expression changes to patterns related to EGFR in- that application of this statistic may require similar fil- hibition [41], perturbation of nodes in the EGFR net- tering to retain highly expressed genes. In the case of work [39], and time course dynamics of targeted DNA methylation, these PatternMarker genes also in- therapeutics. In this dataset, CoGAPS analysis of gene clude genes representing driver alterations in resistance. expression data from cetuximab-resistant clones distin- The DNA methylation data did not require filtering guished the patterns for immediate gene expression when applying the PatternMarker statistic since no flat changes and patterns for long-term changes associated pattern was detected. However, transcriptional regula- with acquired resistance. Gene expression signatures for tion by epigenetic alterations or in pathways involves resistance to EGFR inhibitors in two additional cell lines simultaneous co-regulation of multiple genes. This co- (one HNSCC and one NSCLC) from previous studies regulation is reflected in the reuse of genes in CoGAPS [38, 41] were significantly enriched in both types of gene weights associated with each pattern. Therefore, es- CoGAPS patterns. Since these previous resistance signa- timates of pathway dynamics from transcriptional data tures were learned from case-control studies without require accounting for all genes with gene set enrich- multiple time point measurements, we concluded that ment statistics instead of the PatternMarker statistic. our time course data are instrumental in discriminating Thus, we hypothesize that the PatternMarker statistic is the signatures of immediate therapeutic response from robust for visualization and biomarker identification. On signatures of acquired resistance. the other hand, gene set enrichment of the CoGAPS In spite of the complexities of the data integration, the gene weights corresponding to each pattern and stored weight of each sample in patterns inferred by CoGAPS in the amplitude matrix are essential for characterization reflects the dynamics of the process in each data modal- of functional alterations in pathways. ity. These patterns are learned completely unsupervised from the data and do not require any gene selection or Integrated genomics analysis of time course data comparison between time points relative to any refer- demonstrates that DNA methylation changes follow ence control. The CoGAPS analysis of the time course transcriptional changes, leading to the model where data demonstrates that applying matrix factorization al- methylation stabilizes the resistance phenotype gorithms for genomics can reconstruct signals associated Collecting treated and untreated cells to obtain paired with phenotypes from time course, omics data. The measurements of methylation and gene expression en- genes associated with CoGAPS patterns had weights that abled us to evaluate whether changes in DNA methyla- were non-zero in multiple patterns. The PatternMarker tion impact gene expression. Including a PBS control at statistic [26] enable further selection of the genes that every time point also enabled the discrimination of the are uniquely associated with each pattern. Creating a changes that result from an adaptive response to therapy heatmap of the genomics profiles for these genes en- from changes that result from maintaining cells in cul- abled novel, heatmap-based visualization of the temporal ture. CoGAPS analysis of DNA methylation data denotes dynamics in the omics data. CoGAPS analysis of gene only changes associated with acquired resistance, in con- expression data contains a flat pattern (EP5), which in- trast to the immediate expression changes observed with cludes all highly expressed genes. These genes may also cetuximab treatment. Thus, while therapeutic response change in association with the experimental conditions, can drive massive changes in gene expression, only the Stein-O’Brien et al. Genome Medicine (2018) 10:37 Page 17 of 22 subset of expression changes associated with the devel- gene regulation, associated with acquired therapeutic resist- opment of resistance have corresponding epigenetic sig- ance. Future investigation into the chromatin remodeling natures, suggesting that the methylation landscape is mechanisms will also test whether chromatin alterations important for the development of acquired resistance. follow the changes in expression and occur in combination The CoGAPS patterns in gene expression that are asso- with altered methylation patterns to drive epigenetic regula- ciated with acquired cetuximab resistance gradually tion of resistance. change over the time course (EP1 and EP2). On the other hand, the CoGAPS patterns for DNA methylation changes have a sharp transition at the generation at Time course data are essential to distinguish clonal which resistance is acquired (CTX-G4). These patterns outgrowth from transcriptional rewiring that give rise to (MP2 and MP3) reflect a delayed but more rapid change stable acquired resistance to therapy in DNA methylation. The time delays between alter- Besides the immediate changes in gene expression ations in DNA methylation and gene expression pose a followed by the gradual methylation switch, it is also in- further computational challenge for integrated, time teresting to note these effects in the proliferation rates course genomics analyses. The vast majority of inte- during the 11 weeks of treatment. Initially, the prolifera- grated analysis algorithms assume one-to-one mapping tion of the population of cetuximab-sensitive cells is of genes in different data platforms or seek common pat- slower when compared to the untreated controls, reflect- terns or latent variables across them [53]. Such ap- ing therapy effectiveness. Early and progressively, the proaches would fail to capture the early changes from cells develop molecular changes to overcome the EGFR cetuximab treatment that impact only gene expression; blockade. However, this process starts in just a small time delays between DNA methylation and gene expres- number of clones and the increase in proliferation is still sion patterns and different gene usage in each pattern. It not enough to surpass the growth rate of the untreated is essential to develop new integrated algorithms to sim- cells. As soon as the population of resistant cells is larger ultaneously distinguish both patterns that are shared than the number of sensitive cells, the proliferation rate across data types and that are unique to each platform. is now higher than in the untreated controls. At some For time course data analysis, these algorithms must also point, we observe the stabilization of the proliferation model regulatory relationships that may give rise to tim- rates in the cetuximab-treated cells, probably due to the ing delays, such as epigenetic silencing of gene expres- fact the culture is now dominated by the population of sion. However, as we observed with the unanticipated resistant clones. Although stable, the proliferation rate changes in DNA methylation following and not preced- of these resistant clones is significantly higher than that ing gene expression, they must also consider delays of the untreated cells. This increased proliferation rate is resulting from larger phenotypic changes such as the consistent with the rapid increase in tumor volume ob- stability of the therapeutic resistant phenotype. served clinically once patients develop resistance to ther- The relative timing of change in DNA methylation and apy. Tracing the increase in the population of resistant gene expression is consistent with previous observations cells and their proliferation rates in vivo requires novel that gene expression changes precede DNA methylation al- techniques to biopsy or image tumors at intermediate INK4A terations in genes critical for cancer progression. P16 time points of treatment. and GSTP1 are tumor suppressor genes for which In a recent study, gene expression changes were asso- transcription silencing was found to occur prior to DNA ciated with a transient resistant phenotype present in hypermethylation and chromatin changes. The temporal melanoma cell lines before vemurafenib administration delay observed between expression and methylation [56]. Once the melanoma cells were exposed to the drug, patterns in our time course provides transcriptome-wide additional changes in gene expression are detected and evidence of this phenomenon. Specifically, that epigenetic are later accompanied by changes in chromatin structure changes are necessary to stabilize gene expression aberrant [56]. These findings, together with our time course ob- profile and will be followed by modification into a silenced servations, suggest that in the heterogeneous tumor en- methylation state, resulting in tumor progression [54, 55]. vironment the existence of some cells expressing specific Our integrated RNA-seq and DNA methylation analysis marker genes can trigger cellular reprogramming as corroborates the fact that gene expression changes occur soon as the targeted therapy is initiated. Upon drug ad- earlier to epigenetic alterations and suggests that DNA ministration, the number of genes with aberrant expres- methylation is essential to maintain the changes in gene ex- sion increases and is followed by other epigenetic and pression in this acquired cetuximab resistance model. genetic changes that will shift the transient resistant Additional time course data tracing other in vitro and in state into a stable phenotype. This finding on acquired vivo models of HNSCC are essential to generalize the rela- resistance development could dramatically change the tive timing of molecular changes, and thus mechanisms of course of treatment with targeted therapeutic agents. Stein-O’Brien et al. Genome Medicine (2018) 10:37 Page 18 of 22 The precise characterization of resistant gene signatures inhibition resulted in phospho-ERK decreased expres- and their timing are crucial to determine the correct sion that was restored when FGF2 was added to the point during the patients’ clinical evolution to introduce culture, suggesting that an autocrine FGF-FGFR path- alternative therapeutic strategies. This way, secondary way is one of the mechanisms of resistance to gefitinib. interventions would start before the stable resistant However, the cell lines evaluated in both studies were phenotype is spread among the tumor cells resulting in intrinsically resistant to gefitinib. In our model, we in- prolonged disease control and substantial increase in duced resistance to cetuximab and observed FGFR1 overall survival. gain of expression and significant anti-correlation with the DNA methylation. We additionally evaluated the expression of other FGFRs and FGFs that were identi- DNA methylation changes in FGFR1 are associated with fied by the PatternMarker statistic. Although it is not changes in FGFR signaling found as a PatternMarker in our analysis pipeline, Among the genes we identified with the canonical rela- FGF2 is upregulated in the cetuximab generations when tionship between expression and methylation, FGFR1 compared to the PBS generations as observed with present with increased gene expression accompanied by FGFR1 (Additional file 1: Figure S14). Thus, our data loss of CpG methylation. FGFR1 is a receptor tyrosine corroboratesand extendsthispreviousevidencefrom kinase that regulates downstream pathways, such as intrinsically resistant lines that one of the mechanisms PI3K/AKT, and RAS/MAPK, which are also regulated driving resistance to EGFR inhibition is the FGF-FGFR by EGFR [57]. Its overexpression has been previously autocrine pathway. This observation adds another evi- associated with resistance to EGFR inhibitors in other dence that the computational approach used in this cancer types including HNSCC [27–29]. To our know- study is robust once it is capable of identifying mecha- ledge this is the first study showing that epigenetic al- nisms previously described in other models resistant to terations are associated with changes in FGFR1 EGFR blockade. expression in HNSCC during the development acquired Our previously developed bioinformatics algorithms cetuximab resistance. FGFR1 upregulation combined for the identification of gene expression and epigenetic with promoter hypomethylation was previously de- patterns progression over time proved to be consistent, scribed in rhabdomyosarcomas [58]. Other studies de- since they also detected canonical changes found to be scribed that FGFR1 increased levels is a common driving this mechanism among innumerous new poten- feature in different tumor types, such as glioblastoma tial candidates for acquired resistance. The integrated [59] and cancers of the breast [60], lung [61], prostate computational analysis was possible due to an experi- [62], bladder [63], ovarian [64], colorectal [27], and mental approach developed to account for molecular HNSCC [29, 65, 66]. FGFR1 is involved in resistance changes due to adaptive responses to the culturing sys- mechanisms against EGFR inhibitors [27, 46–48], such tem and the immediate addition of cetuximab. Here, we as cetuximab and gefitinib. Together, the TCGA and present a novel integrated analysis protocol to evaluate Bossi et al. dataset analyses corroborate our findings molecular changes measured by different high through- that FGFR1 gene expression is regulated by epigenetic put techniques over a prolonged time of treatment with changes in HNSCC. Altogether, the epigenetic alter- an FDA-approved targeted therapeutic agent. The lack ation of FGFR1 represents a candidate biomarker of re- of in vivo experimentation to validate our findings was sistance to cetuximab and further studies are critical to compensated by the analysis of two public datasets of identify combination therapies for HNSCC patients that HNSCC, showing that our in vitro findings were also develop acquired cetuximab resistance. present in patients’ samples. Our findings, together with The increased levels of FGFRs and FGFs are believed Marshall et al. [47] and Marek et al. [46], are a strong to play a role in an autocrine mechanism in HNSCC evidence that FGFR1 plays a crucial role in a significant and NSCLC cell lines with intrinsic resistance to the proportion of cases that are resistant to cetuximab or EGFR inhibitor, gefitinib. Using publicly available gene gefitinib. The translational implications are notable expression microarray datasets, Marshall et al. [47]and since FGFR1 inhibition can be used in combination Marek et al. [46] verified concomitant increased levels with EGFR blockade to retard acquired resistance or of FGFRs and their specific FGFs ligands. Particularly, overcome intrinsic resistance. It is important to men- FGFR1 and FGF2 upregulation was observed in the tion that FGFR inhibitors are being currently evaluated same resistant cell lines and hypothesized to be the by clinical trials and could soon become a potential mechanism behind resistance. This was corroborated new therapeutic option for many cancer patients [57]. by functional experiments showing that cells treated Future work evaluating how these combinations impact with pan-FGFR inhibitor were less prone to anchorage- thetimingofacquiredresistance are essential to independent growth. Also, FGF2 silencing or FGFR1 determine the molecular mechanisms that shift Stein-O’Brien et al. Genome Medicine (2018) 10:37 Page 19 of 22 dominant signaling pathways in cancer and thereby when treated with cisplatin [68]. These data and the in- drive resistance. trinsic sensitivity of SCC25 to EGFR inhibition suggests that therapeutic resistance results from random selection of a pre-existing resistant clone. The heterogeneity in CoGAPS gene signatures from the SCC25 time course methylation profiles reflected the complexity of the resist- model are associated with molecular changes in ance mechanisms that can arise from combination therap- additional in vitro models and human tumor data ies in heterogeneous tumors. Future work extending these The main limitation of the current study was the use of protocols to in vivo models is essential to determine the a single cell line model. SCC25 is intrinsically sensitive role of the microenvironment in inducing therapeutic to cetuximab and from this single cell line model, we resistance. Developing in vivo models with acquired thera- generated two groups of samples (CTX and PBS genera- peutic resistance presents numerous technical challenges tions) over the course of 11 weeks. High-throughput that must first be addressed before such time course measurements and analysis were performed for a total of sampling is possible [9]. Pinpointing precise molecular 22 samples. The collection of multiple data points in the predictors of therapeutic resistance will facilitate the analysis had to be accounted for when determining the identification of unprecedented biomarkers and reveal the number of cell lines to be included in the study. We mechanisms by which to overcome acquired therapeutic nonetheless compared our data to gene signatures from resistance to most therapies used to treat cancer. the other isogenic HNSCC resistant model 1CC8 [10], an independent resistant model to an EGFR inhibitor in NSCLC [38], and human tumor data from HNSCC pa- Conclusions tients before cetuximab treatment [42]. Besides the By developing a novel bioinformatics pipeline for inte- number of samples, we also had to take into consider- grated time course analysis, we measured the changes ation the potential batch and technical effects of broad in gene expression and DNA methylation during the cross-platform profiling. Nevertheless, the analysis of progression from an intrinsic cetuximab responsive pretreatment HNSCC patient samples from TCGA [36] state to the acquired resistant phenotype using an in and another study [42] confirmed that our finding that vitro HNSCC cell line model. Specifically, this pipeline FGFR1 is upregulated and demethylated in HNSCC and includes: (1) CoGAPS analysis of each platform inde- associated with acquired resistance to cetuximab is also pendently; (2) gene selection with the PatternMarker a mechanism involved in intrinsic resistance to the tar- statistic for visualization and CoGAPS gene set analysis geted therapy. of the CoGAPS gene profiles for pathway analysis; (3) comparisons of patterns to known phenotypes infer Transcriptional heterogeneity is critical in acquired their relative timing; (4) anti-correlation between DNA therapeutic resistance, and future time course single cell methylation patterns and gene expression to infer epi- data will pinpoint precise molecular predictors of genetically regulated genes; and (5) evaluation of Pat- therapeutic resistance ternMarker genes and projection of the CoGAPS gene The in vitro protocol for time course sampling developed profiles to learn relevance of inferred gene signatures in in this study has the additional advantage of aggregating new datasets. This pipeline revealed massive changes in potentially heterogeneous mechanisms of resistance in- gene expression and identified and discriminated the creasing the signal of changes in any cetuximab-resistant different patterns associated with resistance or cell cul- sub-clone. For example, we observed demethylation and turing conditions. This analysis demonstrates that com- overexpression of FGFR1 in the pooled cells, but only a pressed sensing matrix factorization algorithms can single stable clone generated from the same SCC25 cell identify gene signatures associated with the dynamics line in a previous study (CTXR10) had upregulation of of phenotypic changes from genomics data collected FGFR1 [31]. This finding suggests that tumor heterogen- over the time course. In this case, the gene expression eity also plays a role in acquired resistance to targeted patterns relevant to resistance were later followed by therapies and enables different pathways to be used to epigenetic alterations. Our main conclusion is that bypass the silenced target within the same tumor. using our bioinformatics approach we are able to deter- Heterogeneity of SCC25 cetuximab-resistant clones has mine that the resistant phenotype is driven by gene ex- been observed previously [31]. Recent single cell RNA-seq pression changes that would confer the cancer cells data of SCC25 has shown that there is considerable tran- adaptive advantages to the treatment with cetuximab. scriptional heterogeneity in this cell line before treatment Finally, the integrated analysis show that the stability of [67]. Other cancer therapies are influenced by heterogen- the resistant state is dependent on epigenetic changes eity and outgrowth of resistant clones, as was observed in that will make these new gene signatures heritable to single cell clones isolated from the HNSCC cell line FaDu expand the phenotype to the daughter cells. The Stein-O’Brien et al. Genome Medicine (2018) 10:37 Page 20 of 22 bioinformatics pipeline we developed is also significant acid; RNA-seq: Ribonucleic acid sequencing; rRNA: Ribosomal ribonucleic acid; SPFS: Short-progression-free survival; STR: Short tandem repeat; to clinical practice, since it pointed the time course of TCGA: The Cancer Genome Atlas; TFAP2A: Transcription factor AP-2 alpha molecular changes associated with acquired cetuximab resistance and suggests that the resistant phenotype Acknowledgements can be reversed if alternative interventions are intro- We thank JHMI Deep Sequencing & Microarray Core and SKCCC Microarray duced before epigenetic alterations to the genes driving Core Facility on performing and providing advice on RNA-seq and DNA methylation hybridization arrays, respectively; S. Boca, B. Kerr, S. Floor, C. Mak, acquired resistance. Moreover, the computational T. Ou, D. Sidransky, L. M. Weiner, F. Zamuner, K. Zambo, and members of approach we describe here can be applied to time NewPISlack for critical comments and feedback during the preparation of course studies using other tumor type models and the manuscript. In memory of Elizabeth Klein. targeted therapies. Funding This work was supported by NIH Grants R01CA177669, R21DE025398, Additional files P30CA006973, R01DE017982, SPORE P50DE019032, and the Johns Hopkins University Catalyst Award. Additional file 1: Supplemental Figures S1 - S14. Figure S1 - Time course approach to induce resistance to cetuximab and measure gene Availability of data and materials expression and DNA methylation changes; Figure S2 - Anchorage- Unless otherwise specified, all genomics analyses were performed in R and independent growth of cetuximab generation 10 (CTX-G10); Figure S3 - code for these analyses is available from https://sourceforge.net/projects/ Heatmap and hierarchical clustering of gene expression values in 11 gen- scc25timecourse. All transcriptional and epigenetic data of the cell lines from erations of SCC25 cells treated with PBS as control (black columns) and this paper have been deposited in GEO (GSE98815). DNA methylation data with 100 nM of cetuximab (red columns) to acquire resistance.; Figure of head and neck tumors after treatment with cetuximab are available in S4 - Time course gene expression comparison to previously known gene GEO (GSE110996). signatures of resistance to EGFR inhibitors; Figure S5 - Expected gene ex- pression values for genes in each CoGAPS pattern inferred from gene ex- pression data over generations of PBS control (black lines) or treatment Authors’ contributions with 100 nM of cetuximab (red lines); Figure S6 - Heatmap of gene set GS, LTK, SL, CHC, and EJF planned, designed and wrote the manuscript with analysis scores for targets of transcription factors in theEGFR network, tar- input from all authors. GS, LTK, SL, MT, CHC, and EJF contributed to the gets of the AP-2alpha transcription factors associated with cetuximabre- development of methodology. SS provided DNA from HNSCC patients sponse, and cetuximab resistance signatures in CoGAPS patterns; Figure following cetuximab treatment for DNA methylation analysis. GS, SL, and EJF S7 - Heatmaps of Pearson correlation coefficients between CoGAPS gene performed analysis and interpretation of data (e.g. computational analysis). expression and DNA methylation patterns; Figure S8 - Gene expression RR, HO, HC, MC, AF, LVD, SS, JA, and DAG participated in development of heatmap for the time course experiment vs. single cell resistant clones methodology and provided technical and material support. RR, LVD, EI, and experiment; Figure S9 - DNA methylation heatmap for the time course DAG participated in the review and/or revision of the manuscript. All authors experiment vs. single cell resistant clones experiment; Figure S10 - Mi- discussed the data and contributed to the manuscript preparation. CHC and croscopy images of cetuximab single cell resistant clones (CTXR); Figure EJF instigated and supervised the project. All authors read and approved the S11 - Epigenetically regulated pattern marker genes associated with re- final manuscript. sistance presenting significant anti-correlation between gene expression and DNA methylation in thecetuximab single cell resistant clones (CTXR); Figure S12 - Proliferation assay of cetuximab resistant single cell clones Ethics approval and consent to participate (CTXR); Figure S13 - Bar plot of FGFR1 expression in human tumor sam- HNSCC samples and patient information collection were approved by the ples pre (black) and post (red) treatment with cetuximab; Figure S14 Independent Ethics Committee and the Belgian Health Authorities (ID: 2008/ - Heatmap of FGF family genes methylation (A) and expression (B) in time 20MAR/090) and were conducted in accordance with the Declaration of coursedata for treated and control samples. (PDF 22593 kb) Helsinki (October 2000). It was prospectively planned to perform translational research and patients gave their informed consent for repeated biopsies. Additional file 2: Table S1. PatternMarker genes for the expression and methylation CoGAPS patterns and list of genes with significant anti- correlation between gene expression and methylation. (XLSX 117 kb) Competing interests Additional file 3: Table S2. P values of MSigDB hallmark pathways for The authors declare that they have no competing interests. CoGAPS expression patterns calculated with the permutation-based statistic in calcCoGAPSStat. (XLSX 99 kb) Additional file 4: Table S3. P values of MSigDB hallmark pathways Publisher’sNote for CoGAPS DNA methylation patterns calculated with the permutation- Springer Nature remains neutral with regard to jurisdictional claims in published based statistic in calcCoGAPSStat. (XLSX 79 kb) maps and institutional affiliations. Author details Abbreviations Institute of Genetic Medicine, Johns Hopkins University, Baltimore, MD, USA. AKT: AKT serine/threonine kinase; ATCC: American Type Culture Collection; Department of Oncology, Sidney Kimmel Comprehensive Cancer Center, CoGAPS: Coordinated Gene Activity in Pattern Sets; CTX: Cetuximab; Johns Hopkins University, Baltimore, MD, USA. Department of Head and CTXR: Single cell cetuximab resistant clone; DNA: Deoxyribonucleic acid; Neck-Endocrine Oncology, Moffitt Cancer Center, Tampa, FL, USA. EGFR: Endogenous growth factor receptor; FDA: Food and Drug Department of Otorhinolaryngology-Head and Neck Surgery, Keio University Administration; FGFR1: Fibroblast growth factor receptor 1; GEO: Gene School of Medicine, Tokyo, Japan. Department of Surgery - Otolaryngology– Expression Omnibus; GSTP1: Glutathione S-Transferase pi 1; Head and Neck Surgery, University of Utah, |Salt Lake City, UT, USA. Head GWCoGAPS: Genome Wide Coordinated Gene Activity in Pattern Sets; and Neck Surgery Unit, St Luc University Hospital, Brussels, Belgium. HNSCC: Head and neck squamous cell carcinoma; HRAS: HRAS proto- Laboratory of Systems Biology and Computational Genetics, Vavilov Institute oncogene; JHMI: Johns Hopkins Medical Institutions; LPFS: Long-progression- of General Genetics, Russian Academy of Sciences, Moscow, Russia. free survival; MAPK: Mitogen activated kinase-like protein; NSCLC: Non-small- Department of Surgery, UC San Diego Moores Cancer Center, La Jolla, CA, cell lung cancer; PBS: Phosphate buffered saline; PI3K: Phosphoinositide-3- USA. Department of Otolaryngology-Head and Neck Surgery, Johns Hopkins kinase regulatory subunit 1; RIN: RNA integrity number; RNA: Ribonucleic University, Baltimore, MD, USA. Stein-O’Brien et al. Genome Medicine (2018) 10:37 Page 21 of 22 Received: 23 January 2018 Accepted: 1 May 2018 24. Ochs MF, Fertig EJ. Matrix factorization for transcriptional regulatory network inference. IEEE Symp Comput Intell Bioinforma Comput Biol Proc. 2012;2012:387–96. 25. Fertig EJ, Markovic A, Danilova LV, Gaykalova DA, Cope L, Chung CH, et al. Preferential activation of the hedgehog pathwaybyepigenetic modulations in References HPV negative HNSCC identified with meta-pathway analysis. PLoS One. 2013; 1. Sawyers C. Targeted cancer therapy. Nature. 2004;432:294–7. e78127:8. 2. Hyman DM, Taylor BS, Baselga J. Implementing Genome-Driven Oncology. 26. Stein-O’Brien GL, Carey JL, Lee WS, Considine M, Favorov AV, Flam E, et al. Cell. 2017;168:584–99. PatternMarkers & GWCoGAPS for novel data-driven biomarkers via whole 3. Engelman JA, Zejnullahu K, Mitsudomi T, Song Y, Hyland C, Park JO, et al. transcriptome NMF. Bioinformatics. 2017;33:1892–4. MET amplification leads to gefitinib resistance in lung cancer by activating 27. Azuma K, Kawahara A, Sonoda K, Nakashima K, Tashiro K, Watari K, et al. FGFR1 ERBB3 signaling. Science. 2007;316:1039–43. activation is an escape mechanism in human lung cancer cells resistant to 4. Hata AN, Niederst MJ, Archibald HL, Gomez-Caraballo M, Siddiqui FM, afatinib, a pan-EGFR family kinase inhibitor. Oncotarget. 2014;5:5908–19. Mulvey HE, et al. Tumor cells can follow distinct evolutionary paths to 28. Bertotti A, Papp E, Jones S, Adleff V, Anagnostou V, Lupo B, et al. The become resistant to epidermal growth factor receptor inhibition. Nat Med. genomic landscape of response to EGFR blockade in colorectal cancer. 2016;22:262–9. Nature. 2015;526:263–7. 5. Misale S, Yaeger R, Hobor S, Scala E, Janakiraman M, Liska D, et al. 29. Koole K, Brunen D, van Kempen PMW, Noorlag R, de Bree R, Emergence of KRAS mutations and acquired resistance to anti-EGFR therapy Lieftink C, et al. FGFR1 is a potential prognostic biomarker and therapeutic in colorectal cancer. Nature. 2012;486:532–6. target in head and neck squamous cell carcinoma. Clin Cancer Res. 6. Pietrantonio F, Vernieri C, Siravegna G, Mennitto A, Berenato R, Perrone F, et al. 2016;22:3884–93. Heterogeneity of acquired resistance to anti-EGFR monoclonal antibodies in 30. Hatakeyama H, Cheng H, Wirth P, Counsell A, Marcrom SR, Wood CB, et al. patients with metastatic colorectal cancer. Clin Cancer Res. 2017;23:2414–22. Regulation of heparin-binding EGF-like growth factor by miR-212 and 7. Vincenzi B, Zoccoli A, Pantano F, Venditti O, Galluzzo S. Cetuximab: from acquired cetuximab-resistance in head and neck squamous cell carcinoma. bench to bedside. Curr Cancer Drug Targets. 2010;10:80–95. PLoS One. 2010;e12702:5. 8. Boeckx C, Weyn C, Vanden Bempt I, Deschoolmeester V, Wouters A, 31. Cheng H, Fertig EJ, Ozawa H, Hatakeyama H, Howard JD, Perez J, et al. Specenier P, et al. Mutation analysis of genes in the EGFR pathway in Head Decreased SMAD4 expression is associated with induction of epithelial-to- and Neck cancer patients: implications for anti-EGFR treatment response. mesenchymal transition and cetuximab resistance in head and neck BMC Res Notes. 2014;7:337. squamous cell carcinoma. Cancer Biol Ther. 2015;16:1252–8. 9. Quesnelle KM, Wheeler SE, Ratay MK, Grandis JR. Preclinical modeling of 32. Wang K, Singh D, Zeng Z, Coleman SJ, Huang Y, Savich GL, et al. MapSplice: EGFR inhibitor resistance in head and neck cancer. Cancer Biol Ther. 2012; Accurate mapping of RNA-seq reads for splice junction discovery. Nucleic 13:935–45. Acids Res. 2010;38:e178. 10. Wheeler DL, Huang S, Kruser TJ, Nechrebecki MM, Armstrong EA, Benavente 33. Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data S, et al. Mechanisms of acquired resistance to cetuximab: role of HER (ErbB) with or without a reference genome. BMC Bioinformatics. 2011;12:323. family members. Oncogene. 2008;27:3944–56. 34. Fortin J-P, Labbe A, Lemire M, Zanke BW, Hudson TJ, Fertig EJ, et al. 11. Narayan M, Wilken JA, Harris LN, Baron AT, Kimbler KD, Maihle NJ. Functional normalization of 450k methylation array data improves Trastuzumab-induced HER reprogramming in “resistant” breast carcinoma replication in large cancer studies. Genome Biol. 2014;15:503. cells. Cancer Res. 2009;69:2191–4. 35. Bidaut G. Interpreting and comparing clustering experiments through graph 12. Smyth GK. Linear models and empirical Bayes methods for assessing visualization and ontology statistical enrichment with the ClutrFree package. differential expression in microarray experiments. Stat Appl Genet Mol Biol. In: Ochs MF, Casagrande JT, Davuluri RV, editors. Biomed. Inform. Cancer 2004;3:Article3. Res. Boston, MA: Springer US; 2010. p. 315–33. http://link.springer.com/10. 13. Ernst J, Nau GJ, Bar-Joseph Z. Clustering short time series gene expression 1007/978-1-4419-5714-6_19. Accessed 23 Jan 2018. data. Bioinforma Oxf Engl. 2005;21(Suppl 1):i159–68. 14. Lin D, Shkedy Z, Yekutieli D, Burzykowski T, Göhlmann HWH, De Bondt A, et 36. Lawrence MS, Sougnez C, Lichtenstein L, Cibulskis K, Lander E, Gabriel SB, al. Testing for trends in dose-response microarray experiments: a et al. Comprehensive genomic characterization of head and neck squamous comparison of several testing procedures, multiplicity and resampling-based cell carcinomas. Nature. 2015;517:576–82. inference. Stat Appl Genet Mol Biol. 2007;6:Article26. 37. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: a knowledge-based approach for 15. Aryee MJ, Gutiérrez-Pabello JA, Kramnik I, Maiti T, Quackenbush J. An interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. improved empirical bayes approach to estimating differential gene 2005;102:15545–50. expression in microarray time-course data: BETR (Bayesian Estimation of 38. Coldren CD. Baseline gene expression predicts sensitivity to gefitinib in Temporal Regulation). BMC Bioinformatics. 2009;10:409. non-small cell lung cancer cell lines. Mol Cancer Res. 2006;4:521–8. 16. Liao JC, Boscolo R, Yang Y-L, Tran LM, Sabatti C, Roychowdhury VP. Network 39. Fertig EJ, Ren Q, Cheng H, Hatakeyama H, Dicker AP, Rodeck U, et al. Gene component analysis: reconstruction of regulatory signals in biological expression signatures modulated by epidermal growth factor receptor systems. Proc Natl Acad Sci U S A. 2003;100:15522–7. activation and their relationship to cetuximab resistance in head and neck 17. Bonneau R, Reiss DJ, Shannon P, Facciotti M, Hood L, Baliga NS, et al. The squamous cell carcinoma. BMC Genomics. 2012;13:160. Inferelator: an algorithm for learning parsimonious regulatory networks from systems-biology data sets de novo. Genome Biol. 2006;7:R36. 40. McCall MN, Bolstad BM, Irizarry RA. Frozen robust multiarray analysis (fRMA). 18. Ernst J, Vainas O, Harbison CT, Simon I, Bar-Joseph Z. Reconstructing Biostat Oxf Engl. 2010;11:242–53. dynamic regulatory maps. Mol Syst Biol. 2007;3:74. 41. Fertig EJ, Ozawa H, Thakar M, Howard JD, Kagohara LT, Krigsfeld G, et al. 19. Seok J, Xiao W, Moldawer LL, Davis RW, Covert MW. A dynamic network of CoGAPS matrix factorization algorithm identifies transcriptional changes in transcription in LPS-treated human subjects. BMC Syst Biol. 2009;3:78. AP-2alpha target genes in feedback from therapeutic inhibition of the EGFR 20. Naegle KM, Welsch RE, Yaffe MB, White FM, Lauffenburger DA. MCAM: multiple network. Oncotarget. 2016;7:73845–64. clustering analysis methodology for deriving hypotheses and insights from 42. Bossi P, Bergamini C, Siano M, Cossu Rocca M, Sponghini AP, Favales F, high-throughput proteomic datasets. PLoS Comput Biol. 2011;7:e1002119. et al. Functional genomics uncover the biology behind the responsiveness of head and neck squamous cell cancer patients to cetuximab. Clin Cancer 21. Fernández MA, Rueda C, Peddada SD. Identification of a core set of Res. 2016;22:3961–70. signature cell cycle genes whose relative order of time to peak expression is conserved across species. Nucleic Acids Res. 2012;40:2823–32. 43. Schmitz S, Bindea G, Albu RI, Mlecnik B, Machiels J-P. Cetuximab 22. Wise A, Bar-Joseph Z. SMARTS: reconstructing disease response networks promotes epithelial to mesenchymal transition and cancer associated from multiple individuals using time series gene expression data. fibroblasts in patients with head and neck cancer. Oncotarget. Bioinforma Oxf Engl. 2015;31:1250–7. 2015;6:34288–99. 23. Fertig EJ, Ding J, Favorov AV, Parmigiani G, Ochs MF. CoGAPS: an R/C++ 44. Fortin J-P, Triche TJ, Hansen KD. Preprocessing, normalization and package to identify patterns and biological process activity in transcriptomic integration of the Illumina HumanMethylationEPIC array with minfi. data. Bioinforma Oxf Engl. 2010;26:2792–3. Bioinformatics. 2017;33:558–60. Stein-O’Brien et al. Genome Medicine (2018) 10:37 Page 22 of 22 45. Mootha VK, Lindgren CM, Eriksson K-F, Subramanian A, Sihag S, Lehar J, 67. Puram SV, Tirosh I, Parikh AS, Patel AP, Yizhak K, Gillespie S, et al. Single-cell et al. PGC-1α-responsive genes involved in oxidative phosphorylation transcriptomic analysis of primary and metastatic tumor ecosystems in head are coordinately downregulated in human diabetes. Nat Genet. 2003; and neck cancer. Cell. 2017;171:1611–1624.e24. 34:267–73. 68. Niehr F, Eder T, Pilz T, Konschak R, Treue D, Klauschen F, et al. Multilayered 46. Marek L, Ware KE, Fritzsche A, Hercule P, Helton WR, Smith JE, et al. Fibroblast omics-based analysis of a head and neck cancer model of cisplatin growth factor (FGF) and FGF receptor-mediated autocrine signaling in non- resistance reveals intratumoral heterogeneity and treatment-induced clonal small-cell lung cancer cells. Mol Pharmacol. 2009;75:196–207. selection. Clin Cancer Res. 2018;24:158–68. 47. Marshall ME, Hinz TK, Kono SA, Singleton KR, Bichon B, Ware KE, et al. Fibroblast growth factor receptors are components of autocrine signaling networks in head and neck squamous cell carcinoma cells. Clin Cancer Res. 2011;17:5016–25. 48. Wynes MW, Hinz TK, Gao D, Martini M, Marek LA, Ware KE, et al. FGFR1 mRNA and protein expression, not gene copy number, predict FGFR TKI sensitivity across all lung cancer histologies. Clin Cancer Res. 2014;20:3299–309. 49. Ochs MF, Rink L, Tarn C, Mburu S, Taguchi T, Eisenberg B, et al. Detection of treatment-induced changes in signaling pathways in gastrointestinal stromal tumors using transcriptomic data. Cancer Res. 2009;69:9125–32. 50. Hafner M, Niepel M, Chung M, Sorger PK. Growth rate inhibition metrics correct for confounders in measuring sensitivity to cancer drugs. Nat Methods. 2016;13:521–7. 51. Hill SM, Heiser LM, Cokelaer T, Unger M, Nesser NK, Carlin DE, et al. Inferring causal molecular networks: empirical assessment through a community- based effort. Nat Methods. 2016;13:310–8. 52. Zhu X, Ching T, Pan X, Weissman SM, Garmire L. Detecting heterogeneity in single- cell RNA-Seq data by non-negative matrix factorization. PeerJ. 2017;5:e2888. 53. Tyekucheva S, Marchionni L, Karchin R, Parmigiani G. Integrating diverse genomic data using gene sets. Genome Biol. 2011;12:R105. 54. Bachman KE, Park BH, Rhee I, Rajagopalan H, Herman JG, Baylin SB, et al. Histone modifications and silencing prior to DNA methylation of a tumor suppressor gene. Cancer Cell. 2003;3:89–95. 55. Stirzaker C, Song JZ, Davidson B, Clark SJ. Transcriptional gene silencing promotes DNA hypermethylation through a sequential change in chromatin modifications in cancer cells. Cancer Res. 2004;64:3871–7. 56. Shaffer SM, Dunagin MC, Torborg SR, Torre EA, Emert B, Krepler C, et al. Rare cell variability and drug-induced reprogramming as a mode of cancer drug resistance. Nature. 2017;546:431–5. 57. Babina IS, Turner NC. Advances and challenges in targeting FGFR signalling in cancer. Nat Rev Cancer 2017;17:318–32. 58. Goldstein M, Meller I, Orr-Urtreger A. FGFR1 over-expression in primary rhabdomyosarcoma tumors is associated with hypomethylation of a 5’ CpG island and abnormal expression of the AKT1, NOG, and BMP4 genes. Genes Chromosomes Cancer. 2007;46:1028–38. 59. Rand V, Huang J, Stockwell T, Ferriera S, Buzko O, Levy S, et al. Sequence survey of receptor tyrosine kinases reveals mutations in glioblastomas. Proc Natl Acad Sci U S A. 2005;102:14344–9. 60. Andre F, Bachelot T, Campone M, Dalenc F, Perez-Garcia JM, Hurvitz SA, et al. Targeting FGFR with Dovitinib (TKI258): preclinical and clinical data in breast cancer. Clin Cancer Res. 2013;19:3693–702. 61. Cihoric N, Savic S, Schneider S, Ackermann I, Bichsel-Naef M, Schmid RA, et al. Prognostic role of FGFR1 amplification in early-stage non-small cell lung cancer. Br J Cancer. 2014;110:2914–22. 62. Armstrong K, Ahmad I, Kalna G, Tan SS, Edwards J, Robson CN, et al. Upregulated FGFR1 expression is associated with the transition of hormone- naive to castrate-resistant prostate cancer. Br J Cancer. 2011;105:1362–9. 63. Tomlinson DC, Lamont FR, Shnyder SD, Knowles MA. Fibroblast growth factor receptor 1 promotes proliferation and survival via activation of the mitogen-activated protein kinase pathway in bladder cancer. Cancer Res. 2009;69:4613–20. 64. Gorringe KL, Jacobs S, Thompson ER, Sridhar A, Qiu W, Choong DYH, et al. High-resolution single nucleotide polymorphism array analysis of epithelial ovarian cancer reveals numerous microdeletions and amplifications. Clin Cancer Res. 2007;13:4731–9. 65. Göke F, Bode M, Franzen A, Kirsten R, Goltz D, Göke A, et al. Fibroblast growth factor receptor 1 amplification is a common event in squamous cell carcinoma of the head and neck. Mod Pathol. 2013;26:1298–306. 66. Clauditz TS, Böttcher A, Hanken H, Borgmann K, Sauter G, Wilczak W, et al. Prevalence of fibroblast growth factor receptor 1 (FGFR1) amplification in squamous cell carcinomas of the head and neck. J Cancer Res Clin Oncol. 2018;144:53–61.

Journal

Genome MedicineSpringer Journals

Published: May 23, 2018

References

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off