CDRgator: An Integrative Navigator of Cancer Drug Resistance Gene Signatures

Abstract

Understanding the mechanisms of cancer drug resistance is a critical challenge in cancer therapy. For many cancer drugs, various resistance mechanisms have been identified such as target alteration, alternative signaling pathways, epithelial-mesenchymal transition, and epigenetic modulation. Resistance may arise via multiple mechanisms even for a single drug, making it necessary to investigate multiple independent models for comprehensive understanding and therapeutic application. In particular, we hypothesize that different resistance processes result in distinct gene expression changes. Here, we present a web-based database, CDRgator (Cancer Drug Resistance navigator) for comparative analysis of gene expression signatures of cancer drug resistance. Resistance signatures were extracted from two different types of datasets. First, resistance signatures were extracted from transcriptomic profiles of cancer cells or patient samples and their resistance-induced counterparts for >30 cancer drugs. Second, drug resistance group signatures were also extracted from two large-scale drug sensitivity datasets representing ~1,000 cancer cell lines. All the datasets are available for download, and are conveniently accessible based on drug class and cancer type, along with analytic features such as clustering analysis, multidimensional scaling, and pathway analysis. CDRgator allows meta-analysis of independent resistance models for more comprehensive understanding of drug-resistance mechanisms that is difficult to accomplish with individual datasets alone (database URL: http://cdrgator.ewha.ac.kr).

INTRODUCTION

Overcoming drug resistance is one of the most critical challenges in cancer therapy. Cancer drug resistance appears not only in patients who have genetic factors interfering with drug actions, but is also induced by long-term treatment even in initially responsive patients. Because such acquired resistance leads to significant survival disadvantages, it is important to understand the underlying mechanisms of drug resistance in different patients.

In the early studies on resistance mechanisms, genetic mutations that directly alter drug target proteins (i.e. T790M in EGFR and L1152R in ALK) received significant attention (Choi et al., 2010; Pao et al., 2005). These mutations at gatekeeper residues are a common mechanism of resistance to drugs targeting oncogenic kinases, and many patients develop resistance within a year (Holohan et al., 2013). Recently, mutations other than those in the target protein have also been actively studied including alteration of regulatory regions (Leucci et al., 2018; Melton et al., 2017). Although a number of studies on drug resistance have focused on mutations, it has become clear that various genetic events prior to the acquisition of resistance also play an important role (Leucci et al., 2018).

In addition to mutations in the target, other genetic or epigenetic changes including copy number variation (CNV, Bean et al., 2007), activation of bypass signaling (Niederst and Engelman, 2013), and epigenetic modulations (Hu and Baeg, 2017; Nyce et al., 1993) have been identified as primary causes of resistance. Resistance mechanisms are not limited to alterations in only a few genes, but involve overall cellular changes, e.g. epithelial-mesenchymal–transition (EMT) or transition to cancer stem cells (CSCs), which seem to be important factors (Huang et al., 2015; Phi et al., 2018; Wang et al., 2016b). Because such major changes in cell state and physiology cannot be represented by a few genetic alterations, investigation of genome-scale signatures such as transcriptome profiles is likely to be more informative. Transcriptome profiles reflect changes in both genetic and epigenetic conditions including disease progression, drug perturbation, and drug resistance. Therefore, gene expression signatures are useful markers to investigate drug resistance states, classify distinct resistance mechanisms, and predict treatment outcomes.

Typically, there are two strategies to study genetic features associated with drug sensitivity. One is to use large-scale pharmacogenomic datasets of cancer cell lines as model systems, with the cells being divided into two groups, i.e. sensitive and resistant, and genomic features differentiating the two groups being identified. The other is to culture cancer cells in the presence of a drug, leading to selection of resistant cells. By comparing the molecular profiles of the resistance-induced cells with those of the original ones, the relevant genomic features can be analyzed more explicitly than by using the group-based method because the origins of the cells are the same in the former method.

There are several publicly available pharmacogenomics databases such as Genomics of Drug Sensitivity in Cancer (GDSC, https://www.cancerrxgene.org/)(Yang et al., 2013), Cancer Cell Line Encyclopedia (CCLE, https://portals.broadinstitute.org/ccle)(Barretina et al., 2012), and Cancer Therapeutic Response Portal (CTRP, extended version of CCLE, https://portals.broadinstitute.org/ctrp.v2.1/)(Rees et al., 2016). From studies of these databases, mutations in or expression of individual genes associated with drug sensitivity have been characterized. These features have been found useful in predicting resistance based on the genetic background of patients (Iorio et al., 2016; Seashore-Ludlow et al., 2015). While these datasets are among the largest resources to study drug resistance mechanisms, the sparsity of individual mutations makes it difficult to extract genetic markers related to resistance in a comprehensive manner. Only a small fraction of all acquired resistance events in patients may be explained using the markers extracted from GDSC, CCLE, CTRP, and the combination dataset. For example, though the T790M mutation in the EGFR gene has been well known as one of the most frequent alterations associated with acquired resistance to EGFR-targeted drugs in patients (Ma et al., 2011), this mutation has been found in only two cancer cell lines in the above pharmacogenomics databases. Although the Catalog of Somatic Mutations In Cancers (COSMIC, https://cancer.sanger.ac.uk/cosmic)(Forbes et al., 2017) has some genetic features associated with acquired resistance, there is no systematic database of information regarding acquired resistance yet. In COSMIC, mutation information is available for only 13 genes associated with resistance to 23 drugs. There are other databases unrelated to cancer therapy such as databases of antibiotic or antimicrobial drug resistance such as CARD (Jia et al., 2017) and MEGARes (Lakin et al., 2017), respectively. Notably, cancerDR (Kumar et al., 2013) may be the only database dedicated to cancer drug resistance. However, because the resistance-related features in cancerDR mostly originate from CCLE and COSMIC, it is likely to have limitations similar to those of its original resources.

In the present study, we present a web platform, CDRgator (Cancer Drug Resistance signature navigator), which is among the most comprehensive databases of acquired cancer drug resistance because it amasses both resistanceinduced signatures (S_ind) and group-based signatures (S_grp) of drug sensitivity, with the S_ind defined as the differential gene expression in resistant cells of the same lineage induced due to prolonged culture of a cell line in the presence of drug, and S_grp is defined as that in resistant or sensitive groups of cells from different origins. Due to the limitations in genomic features due to sparse mutations, we focused on extracting gene expression signatures of cancer drug resistance to obtain a comprehensive view of resistance mechanisms. For constructing resistance signatures of resistanceinduced cells (S_ind), we performed extensive manual curations of literature as well as data depositories such as Gene express omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) (Barrett et al., 2013) and ArrayExpress (https://www.ebi.ac. uk/arrayexpress/)(Kolesnikov et al., 2015). Alternatively, we also extracted group-based resistance signatures for cancer drugs from CCLE and CTRP; we refer to these resistance signatures extracted from the two large-scale pharmacogenomics databases as group-based signatures (S_grp). These two types of signatures are complementary to each other, and are expected to provide a more comprehensive view on drug resistance, particularly when there are multiple, independent mechanisms involved in developing resistance to a single drug.

Currently, CDRgator provides 603 resistance signatures for 37 cancer drugs representing more than 26 cancer types in total, and the number of signatures will grow as more data are collected. It allows users to browse resistance signatures based on on drugs or cancer types, and to analyze the similarity between resistance signatures. Additionally, CDRgator has a tool to identify resistance signature-matched gene sets from Kyoto Encyclopedia of Genes and Genomes (KEGG, (Kanehisa, 2004) or Gene ontology (GO) (Carbon et al., 2017) to characterize the biological processes involved in resistance. CDRgator also has the ability to compare the resistance signatures in a database with signatures input by users to filter out drugs expected to be resistant or ineffective. Using an illustrative analysis of EGFR inhibitors, we show the utility of CDRgator in understanding the diverse mechanisms of cancer drug resistance.

MATERIALS AND METHODS

Data sources and processing

To obtain induced drug resistance signatures (Sind), we manually collected datasets containing transcriptomic profiles of resistance-induced cells in GEO and ArrayExpress (Fig. 1). The datasets were filtered based on the following criteria; (1) the presence of matched pairs of sensitive and resistant cells, (2) treatment with monotherapy but not combinational therapy for a specific period of time, (3) gene transcription quantification using RNA-seq or microarray (4) in human cells and not in mouse cells or xenografted mouse cells. The RNA-seq data (.fastq file) were mapped using STAR aligner (Dobin et al., 2013) and quantified using HTSeq (Anders et al., 2015). The list of the collected datasets is available with detailed information in Additional file 1.

Fig. 1. The process of extracting drug resistance signatures. Left, resistance-induced signatures; Right, resistance group signatures.

For identifying the resistant group signature (S_grp), we used drug sensitivity data and gene expression data of cancer cell lines from CTRP and CCLE, respectively. To designate the resistant cell line group for a given drug, we grouped cancer cell lines based on their drug responses. The drug responses were represented as area under curve (AUC) of cell growth at different drug concentrations, and normalized with the maximum area calculated assuming 100% response in given concentration ranges. Then, we performed Z-transformation of logged AUC values (Z-AUC) to specify resistant or sensitive cells. We defined cells as a resistant cell group if their ZAUC was less than -0.5. For all drugs, we generated signatures only in cases with at least three cells in the resistant group.

The expression data were downloaded from GEO or ArrayExpress, and were processed using our protocols for microarray and RNAseq analyses. Background correction and normalization was performed using R. For RNAseq, raw fastq files were downloaded and aligned to the reference genome (GRCh37) using STAR aligner. For RANseq of cancer cell lines, BAM files were downloaded from Genomic Data Commons (GDC). To estimate expression, RNAseq read count was calculated using HTseq. The gene expression signature was obtained using differentially expressed gene (DEG) analysis of the resistant cell group or the resistanceinduced cells. For datasets containing biological replicates, genes with adjusted p-values less than 0.05 were selected as DEGs using the limma package (Ritchie et al., 2015) for microarray, and the DESeq2 package (Love et al., 2014) for RNAseq. For datasets with no biological replicates, we calculated fold change in gene expression values and transformed them to z-score. Genes with absolute fold change in z-score of more than 2 were selected as DEGs.

Data statistics

CDRgator contains a total of 143 S_ind for 37 anti-cancer drugs and 30 cell lines from 16 tissue origins. A set of 499 S_grp in our study represented 19 drugs which were also included in S_ind; in terms of details, 267 S_grp were from 20 tissue-level cancer types (e.g. blood cancer) and 232 S_grp were from 23 disease level-cancer types (e.g. acute lymphoblastic leukemia).

Signature similarity

CDRgator provides similarity analysis which explores how similar two signatures are; an enrichment factor (EF) was developed to measure significance of overlap between the signatures compared to a random occurrence. It represents the number of genes common between two sets divided by the expected overlap. At this time, every gene set is restricted to gene spaces common to respective platforms of a signature pair. EF is calculated using the following equation:

Enrichment factor (A, B) = \(\frac{|A \cap B|+p c}{N \times \frac{|A|}{N} \times \frac{|B|}{N}+p c}=\frac{|A \cap B|+p c}{\frac{|A| \times|B|}{N}+p c}\)

where A and B indicate a set of signature genes, and N represents the number of total genes common to both platforms used; pc is a pseudo count that makes the EF robust at small count levels (5 in this study), similar to additive smoothing. EF has the particular advantage of being robust despite differences in signature size.

Functional annotation analysis

To annotate the biological functions of gene signatures, an over-representation test was conducted using the ‘clusterProfiler’ in the R package (Yu et al., 2012), which implements a hypergeometric test of overrepresentation of a gene set against background. Each signature analysis was performed with pathway terms consisting of 10-500 genes among the KEGG pathways or GO biological processes. As the size of the signature set decreases, the number of significantly enriched pathways decreases. Thus, CDRgator also generates a p-value representing the significance of functional enrichment.

Fig. 2. Web interface design of CDRgator. (A) selection of drug class (left circle) and cancer type (right circle); (B) selection of individual datasets; (C) hierarchical clustering of resistance signatures; (D) 2D mapping of resistance signatures via multidimensional scaling (MDS); (E) pathway analysis of resistance signatures using gene set analysis (GSA).

User interface

The CDRgator web service is freely available at http://cdrgator. ewha.ac.kr. Extraction of signatures from expression data and analysis of CDRgator web platform are implemented in R script (R version 3.4.4). All the generated data were stored in the MySql database (v5.7.22). Control of the data from a database has been implemented in Java (JRE 1.8.0) and we built a user-friendly web interface using JSP and jQuery to conveniently provide the data to the user. In addition, the CDRgator web interface uses html5, css3, SASS, and jQuery to load content intuitively and quickly on a variety of platforms, including desktop and mobile. For rich and intuitive visualization of data, it uses the clustergrammer (Fernandez et al., 2017) which is based on the D3.js JavaScript library to build an advanced web interface with animation. The help documentation provides a video and a brief description of the CDRgator guide.

CDRgator is primarily composed of ‘Browse’ and ‘Analysis’ menus. The ‘Browse’ menu provides the ability to search and browse information regarding resistance signatures. The ‘Analysis’ menu performs signature similarity and functional analysis. Signature similarity analysis calculates how similar resistance signatures are, so that users can identify if signatures of specific cancers, drugs, or their classes show a significant tendency in original spaces or two dimensions through multidimensional scaling. A functional analysis shows a heatmap that is generated from an over-representation test to identify genetic and biological functions of signatures. Finally, significant genes can be visualized using KEGG pathways (Supplementary Fig. S1).

RESULTS

Analysis of EGFR inhibitor resistance-induced signatures

To demonstrate the utility of CDRgator, we performed a case study of induced EGFR inhibitor resistance signatures. There were 22 induced resistance signatures for five different EGFR inhibitors. Because these signatures were generated independent of each other, the following questions could be asked; i) is the resistance pattern dependent on the drug, the cancer type, or both; and ii) do drugs with similar modeof-action (e.g. EGFR inhibitors) share changes in transcriptomic profiles.

Similarity clustering analysis of the 22 S_ind for EGFR inhibitors showed a tight cluster, with mostly erlotinib resistance signatures from head and neck cancer (Fig. 3A). This suggests distinct and reproducible patterns of resistance mechanisms specific to a drug and cancer type, because similar signature patterns were generated from multiple independent experiments. This trend was also confirmed in 2D mapping of the signatures using multidimensional scaling (Fig. 3B). These signatures were distinguishable from those of other cancer types. Overall, the resistance patterns collected in CDRgator appear to be heterogeneous even for the same drug and cancer type (e.g. erlotinib (center green) in nonsmall cell lung cancinoma (NSCLC, green outer ring) in Fig. 3B).

Fig. 3. Analysis of resistance-induced signatures for EGFR inhibitors Illustrative case results available in the ‘Analysis’ menu of CDRgator, generated using resistance-induced signatures for five EGFR inhibitors. Color legends of drug and cancer type are described in Fig. 3C; (A) signature similarity analysis; heat represents all pairwise similarities; hierarchical clustering of similarity is shown at the top and left; (B) multidimensional scaling plot of similarity metrics; outer circles are colored cancer types, and inner circle are drugs; the dashed line indicates a cluster including primarily head and neck cancer signatures; (C) KEGG pathway enrichment analysis; up (red) or down (green) regulated gene sets in individual signatures were tested using each pathway term on the left; heat indicates –log p-value of hypergeometric testing between each pathway gene set and signatures; hierarchical clustering is shown on the top and left.

Pathway analysis showed the heterogeneity in resistance signatures in more detail. In head and neck cancers, biological pathways previously known to be associated with drug resistance were also enriched in our signatures, such as upregulation of NF-kB, NOD-like receptor, TNF, and FoxO signaling pathways. NF-kB signaling is well-known among the survival and resistance pathways against many anticancer drugs including EGFR-targeted drugs (Bentires-Alj et al., 2003; Hertlein and Byrd, 2010; Lagunas and MeléndezZajgla, 2008; Wang et al., 2018). It is known that the NF-kB activating complex drives cancer cells to become resistant against EGFR inhibitors (Blakely et al., 2015), and this signaling pathway could also lead to resistance against EGFR inhibitors designed to target mutated forms of EGFR (Galvani et al., 2015). In colon cancers, activated NF-kB induces drug resistance through the regulation of MDR1 gene expression (Bentires-Alj et al., 2003). EGFR triggers proto-oncogenic signals such as RAS and mitogen activated protein kinase (MAPK) signaling (Wee and Wang, 2017). NOD-like receptors are known to regulate innate immunity by forming inflammasomes, and also by activating MAPK and NF-kB signaling (Saxena and Yeretssian, 2014; Shaw et al., 2010). NOD1/2 downstream signaling was reported to confer drug resistance through RIP2-PAX5 interaction by activating MAPK and NF-kB signaling (Wang et al., 2016a). AKT-PI3K signaling is also known to be activated in EGFR drugresistant cancer cells (Jacobsen et al., 2017; Ma et al., 2016). FOXO is among the downstream genes playing a role in cell survival and apoptosis. Drug sensitivity was restored via inhibition of FOXO in AKT-mediated EGFR inhibitor-resistant lung cancers (Sangodkar et al., 2012). Several of these pathways (i.e. NF-kB, NOD-like receptor, and FOXO signaling) were also frequently upregulated in other resistance signatures (Fig. 3B).

In NSCLC, additional pathways (i.e. PI3K-Akt, ECMreceptor interaction, and focal adhesion) were strongly enriched. These pathways were found to be associated with epithelial–mesenchymal transition (EMT), which has been recently highlighted as major resistance mechanism in both chemotherapy and targeted cancer drug therapy (Chen et al., 2013; Huang et al., 2015; Larue and Bellacosa, 2005; Le Bras et al., 2012; Wang et al., 2016b). As shown in the example with EGFR inhibitors, CDRgator provides rich information on the drug resistance processes in terms of heterogeneity, and dependence on cancer type and drug class.

DISCUSSION

The development of technology coupled with bioinformatics and system biology approaches has enabled the identification of genomic and transcriptome features to predict response and resistance to specific drugs. The wide range of molecular mechanisms involving these features has important implications for understanding and treating resistance. However, inhibiting one pathway can result in a relatively simple escape route for the tumor, and an integrated analysis approach involving various features to overcome or circumvent drug resistance is needed.

CDRgator was developed for studying the acquired resistance mechanisms against cancer drugs. Gene expression signatures were generated from manually collected data regarding cancer drug resistance-induced cells in public datasets. CDRgator allows inspection and comparison of these acquired resistance transcriptomic signatures. Pathway analysis enables functional interpretation of multiple signatures simultaneously. While the datasets available in CDRgator are still sparse and insufficient to provide comprehensive information on drug resistance mechanisms, it may provide a useful clue to investigate drug resistance mechanisms, and to develop the right therapeutic strategy. For example, an analysis of genome-scale transcript signatures provided by CDRgator and recurrent mutations in the regulatory region may provide a variety of data enabling resistance marker identification or mechanistic analysis. The complete datasets in CDRgator are accessible at http://cdrgator.ewha.ac.kr and are also downloadable as mysql database dump.

Note: Supplementary information is available on the Molecules and Cells website (www.molcells.org).