Application of single-cell RNA sequencing in optimizing a combinatorial therapeutic strategy in metastatic renal cell carcinoma
- Kyu-Tae Kim†1,
- Hye Won Lee†2,
- Hae-Ock Lee1, 3,
- Hye Jin Song4,
- Da Eun Jeong5,
- Sang Shin5,
- Hyunho Kim6,
- Yoojin Shin6,
- Do-Hyun Nam5, 7,
- Byong Chang Jeong8,
- David G. Kirsch9, 10,
- Kyeung Min Joo4, 5Email author and
- Woong-Yang Park1, 3, 5Email author
© Kim et al. 2016
Received: 21 December 2015
Accepted: 11 April 2016
Published: 29 April 2016
Intratumoral heterogeneity hampers the success of marker-based anticancer treatment because the targeted therapy may eliminate a specific subpopulation of tumor cells while leaving others unharmed. Accordingly, a rational strategy minimizing survival of the drug-resistant subpopulation is essential to achieve long-term therapeutic efficacy.
Using single-cell RNA sequencing (RNA-seq), we examine the intratumoral heterogeneity of a pair of primary renal cell carcinoma and its lung metastasis. Activation of drug target pathways demonstrates considerable variability between the primary and metastatic sites, as well as among individual cancer cells within each site. Based on the prediction of multiple drug target pathway activation, we derive a combinatorial regimen co-targeting two mutually exclusive pathways for the metastatic cancer cells. This combinatorial strategy shows significant increase in the treatment efficacy over monotherapy in the experimental validation using patient-derived xenograft platforms in vitro and in vivo.
Our findings demonstrate the investigational application of single-cell RNA-seq in the design of an anticancer regimen. The approach may overcome intratumoral heterogeneity which hampers the success of precision medicine.
Clear cell renal cell carcinoma (ccRCC), the most prevalent type of sporadic kidney cancer, is often associated with malignant disease progression and poor therapeutic outcomes . A major underlying genetic alteration in ccRCC is the von Hippel–Lindau (VHL) tumor-suppressor gene, whose deregulation stimulates an oncologic metabolic shift . Signaling pathways involved in this metabolic shift have been proposed as potential therapeutic targets, including epidermal growth factor receptor (EGFR) , vascular endothelial growth factor (VEGFR) , or mammalian target of rapamycin (mTOR) pathways . Although targeting these pathways significantly improved progression-free survival, the outgrowth of drug-resistant clones reduced the clinical efficacy and remains a clinical challenge that must be overcome .
Approximately 30 % of patients with renal cell carcinoma (RCC) are diagnosed with metastases . Metastatic renal cell carcinoma (mRCC) evolves from primary RCC (pRCC) and harbors multiple subpopulations with distinct genomic [7, 8] and transcriptomic [9–11] features. Such intratumoral heterogeneity (ITH) is augmented by spatiotemporal tumor evolution [7, 8, 12]. The existence of disparate admixtures of cancer cells across mRCC and pRCC typically leads to significant differences in their sensitivity to therapies  and rare complete response to targeted molecular agents . Therefore, independent investigation of mRCC and pRCC tumor tissues is essential to comprehend the intratumoral landscape within a patient and ultimately to achieve sustainable therapeutic benefit through a rational drug combination design. The magnitude of ITH in ccRCC has been discerned in detail by recent approaches of sequencing multiregional biopsied specimens [7, 8] or single cells . However, experimental applications of the analyzed ITH signature to effectively eradicate tumor cells have not been extensively investigated.
Accurate identification of disease-causing sequence variants and driving pathways is essential to minimize drug resistance or tumor relapse with targeted therapeutics . Genomic mutations, however, may have limited functional significance as druggable targets, according to previous findings that drug responses can be widely different even in genetically homogeneous cancer cell lines [16, 17]. A systematic assessment from the National Cancer Institute and the Dialogue on Reverse Engineering Assessment and Methods (NCI-DREAM) project has shown that gene expression profiling has the best predictive power among independent profiling platforms, with increasing power upon data integration . In the prediction of drug sensitivity in cancer cells, transcriptome profiling can enhance our understanding of how the cellular mechanism is functionally perturbed in response to drug treatment . Moreover, considering that single targeting agents may eliminate a certain subpopulation of tumor cells while leaving others unharmed, it is necessary to analyze the tumor transcriptome at high resolution to detect drug-resistant clones that may be concealed within ITH. Gene expression profiling at a single-cell resolution may enable modeling of functional heterogeneity and identification of subpopulations with specific drug responsive signatures.
Here, we used single-cell RNA sequencing (scRNA-seq) not only to elucidate transcriptional heterogeneity during the metastatic progression, but also to design an optimized combination of targeted agents against metastatic RCC. On the basis of single-cell transcriptome analysis, identification of cellular subpopulations with distinct activation status of signaling pathways allowed us to postulate a combinational therapeutic regimen with an increased likelihood of covering all potentially targetable cancer cells.
Establishment of patient-derived xenografts from paired pRCC and mRCC
Evolutionary genomic trajectories during tumor progression and metastasis
In recognition of that genomic features were consistently propagated with higher cancer cell fraction (~100 %) through xenograft passaging (Additional file 1: Figure S1), we investigated genomic architectures in the pRCC and mRCC tumors from PDX samples to understand the clonal evolution associated with the spatiotemporal tumor progression. WES analysis of primary and paired metastatic samples revealed that 23.5 % of SSNVs were shared (Additional file 2: Figure S2A). In particular, a VHL D121G mutation was found in both samples with high allele frequencies (~1.0, Additional file 2: Figure S2A and Additional file 3: Table S1), suggesting that it might be a founder event in tumor evolution [7, 8]. Variant allele frequencies (VAF) of the shared SSNVs were typically higher than those of SSNVs exclusively observed in mRCC (38 %) or pRCC (38.6 %) (Additional file 2: Figure S2A). Discordant SSNVs in mRCC and pRCC might result from the gradual increase in point mutations and clonal selection with tumor evolution, as previously reported [7, 8]. In contrast, somatic copy number alterations (SCNAs) in mRCC were similar to those in pRCC (Additional file 2: Figure S2B), with 5q amplifications detected only in pRCC (Additional file 2: Figure S2C). Integrated analyses of WES and aCGH to infer evolutionary trajectories showed that major clones harboring driver mutations were shared at high cellular frequencies, whereas minor subclones were enriched in mRCC (Additional file 4: Figure S3A, B and Additional file 5: Table S2). Overall, the RCC of our patient showed a complex non-linear branching clonal evolution (Additional file 4: Figure S3C) that may become the basis of intratumoral diversity [7, 8, 12]. The genetic complexities might also result in molecular and functional differences between pRCC and mRCC despite their clonal origin, as previously reported [9–11].
Single-cell RNA sequencing and quality assessment for expression profiling
To model the functional heterogeneity and to identify specific subpopulations that are phenotypically relevant to drug responses, we used scRNA-seq to profile single cells from the parental mRCC and PDX mRCC and pRCC (Fig. 1b and see “Methods”). After filtering out poor-quality cells, a total of 116 tumor cells from the parental mRCC (n = 34), PDX-mRCC (n = 36), and PDX-pRCC (n = 46) were used in subsequent analyses (Additional file 6: Figure S4 and Additional file 7: Table S3). When compared to the normal kidney cortex, single cancer cells had much more variable gene expression as shown by the high coefficient of variation for averaged gene expression (Additional file 8: Figure S5A). Nonetheless, housekeeping genes, including glyceraldehyde 3-phosphate dehydrogenase (GAPDH) and beta-actin (ACTB), were stably expressed across single cells with low variation within (Additional file 8: Figure S5A) and across (Additional file 8: Figure S5B) cell populations, suggesting functional significance of the gene expression heterogeneity.
Compared to the normal kidney cortex, parental and PDX cancer cells demonstrated low stromal and high ccRCC gene expression signatures (Additional file 9: Figure S6). Interestingly, the global expression profiles of patient-mRCC and PDX-mRCC single-cell samples were more closely related than those of PDX-pRCC samples (Fig. 1c and Additional file 10: Figure S7). Principal component analysis (PCA) across all samples showed three distinctive main clusters (Fig. 1d); in addition to a cluster of normal kidney cortex, a cluster of parental mRCC cells and PDX samples were separated from a cluster of PDX-pRCC samples. As identified in unsupervised clustering of global expression across samples (Fig. 1c), we found the averaged expression of single cells correlated well with that of their bulk cell population samples (Additional file 10: Figure S7B). Furthermore, multiple regression analysis on the transcriptomes of different sized pools of single cells to those of bulk cell population samples showed better representation of the bulk cell population with increasing number of single cells (Additional file 10: Figure S7C). Similarly, in the PCA plot, bulk cell population samples almost matched single cells (Fig. 1d). As the quality assessments and unsupervised clustering analyses demonstrate the reliability of our data, we further explored the single-cell transcriptomes.
Metastatic and aggressive molecular signatures in mRCC transcriptomes
Consistent with the previous observation that PDX cells reflect their parental tumors with rare contamination of normal stromal cells , we verified discrete global expression profiles of single RCC cells compared with normal kidney expression profiles (Fig. 1c, d, Additional file 9: Figure S6, and Additional file 10: Figure S7). The stromal signature of parental bulk mRCC was higher than in other tumor samples (Additional file 9: Figure S6A), concordant with genomic analysis showing a smaller tumor portion in parental mRCC than PDX samples (Additional file 4: Figure S3B). By comparison, single cells from the parental mRCC all showed low stromal signatures (Additional file 9: Figure S6A) and stably higher ccRCC signatures (Additional file 9: Figure S6B), suggesting selective capture of tumor cells from the bulk cell populations. This might be explained by the experimental approach of using large size (17–25 μm) microfluidic chips to capture single cells based on the histopathologic properties of the RCC tumor cells (Fuhrman Grade-III, average cell size ~20 μm)  to reduce the likelihood of capturing stromal cells.
Prediction of activation of drug target pathways and sensitivity to drug responses
To refine the predictive power of signaling pathway activation for drug sensitivity, we built a ridge regression model  using public gene expression profiles and drug sensitivity data as a training set . For the PDX-mRCC and pRCC bulk cell population samples, we found strong positive correlations between predicted and measured drug sensitivity (Fig. 3c). Drug sensitivity was also estimated in all single cancer cells and in normal kidney cortex as a control; these data suggested afatinib and dasatinib sensitivity of pRCC cells (Fig. 3d), similar to the predicted and measured drug sensitivity of bulk cell populations (Fig. 3b). The prediction of drug sensitivity by the ridge regression model was in agreement with the estimation of drug target pathway activation across single cells (Fig. 3e and Additional file 16: Figure S12). Interestingly, pazopanib and everolimus, which showed no clinical benefit in this patient (Fig. 1a), had minor effects on the viability of PDX-pRCC and PDX-mRCC cells (Fig. 3b and Additional file 17: Figure S13) and signaling pathway activation for their targets mTOR and VEGFR also remained low (Fig. 3a and Additional file 13: Figure S10). Although drug sensitivity of cancer cells was predicted to be higher than that of normal kidney cortex for all of the anticancer agents, the difference in pazopanib sensitivity was small (Additional file 16: Figure S12A). Overall, drug sensitivity predictions showed significant correlations with corresponding signaling pathway activation (Additional file 16: Figure S12B). We estimated the activation status of signaling pathways and drug sensitivity and demonstrated concordance between the predicted signatures and the measured data. These results suggest that molecular targeted therapies can be designed on the basis of prediction signatures obtained from RNA-seq.
Heterogeneity in drug response signatures within a tumor
The heterogeneous drug sensitivity prediction for individual tumor cells suggests the presence of tumor subpopulations with distinct signaling pathway activation and drug sensitivity. Based on the prediction results, we identified four mRCC subpopulations with corresponding signaling pathway activation for EGFR and Src pathways (both active, only EGFR active, only Src active, both inactive; Fig. 4). Only a fraction of cells (23.5 % in parental mRCC; 27.8 % in PDX-mRCC) with activated status for both in EGFR and Src signaling pathways were predicted to be sensitive to both afatinib and dasatinib. A population with inactivated status in both signaling pathways (14.7 % in parental mRCC; 13.9 % in PDX-mRCC) is unlikely to respond to afatinib or dasatinib. The largest proportion of cells (61.8 % in parental mRCC; 58.3 % in PDX-mRCC) had signaling pathway activation in either of the two pathways, suggesting a mutually exclusive response to the corresponding drugs. This classification suggests that co-targeting of EGFR and Src with a combination of afatinib and dasatinib would be an efficient therapeutic strategy without redundancy and with additive antitumor effects (Fig. 4e).
Evaluation of the single-cell analysis-driven therapeutic strategy
Finally, we tested the in vivo performance of the combinatorial drug treatment using the subcutaneous xenograft model of mRCC. As single agents, afatinib (20 mg/kg) or dasatinib (15 mg/kg) had modest inhibitory effects on tumor growth (afatinib, 55 %; dasatinib, 40 %; Fig. 5d). Despite the partial treatment effect, the tumor cells continued to grow and none of the tumor xenografts showed complete response. In combination, afatinib and dasatinib showed a significantly enhanced antitumor effect, inhibiting tumor growth by 78 % (Fig. 5d). Single-agent and combination treatment protocols were well tolerated by the mice, with no weight loss or other signs of acute or delayed toxicity (Fig. 5e). The drug effects were confirmed at the molecular level by complete inhibition of AKT (Fig. 5f) and ERK (Fig. 5g) phosphorylation. Together, co-targeting functionally distinct subpopulations in mRCC PDX cells that were identified by single-cell transcriptome analyses significantly improved treatment outcomes in both in vitro and in vivo investigational models.
The evolution of multiple tumor subclones during tumor progression and metastasis generates intratumoral heterogeneity, which plays a role in intrinsic and acquired treatment resistance to molecular targeted therapies. In mRCC, durable complete responses are rarely achieved with conventional targeted molecular agents . Furthermore, primary tumors and their corresponding metastases frequently show significant differences in therapeutic responses , suggesting that biopsies should be taken upon metastatic relapse for molecular profiling analysis to inform selection of salvage therapies. To date, most molecular profiling has been performed in the primary tumors of RCC patients and therefore failed to recapitulate the metastatic population. To overcome such clinical challenges, therapeutic strategies that efficiently target heterogeneous tumor subpopulations in mRCC must be developed.
In this report, we introduce scRNA-seq as a promising strategy that allows high-resolution analyses of the intratumoral landscape for tumor transcriptomes  and optimization of targeted treatment regimens against metastatic RCC. In our model case, we found that mRCC diverged from pRCC both at genomic and transcriptomic levels, manifesting distinctive genetic aberrations and drug target pathway activation. We combined comprehensive single-cell analyses with high-throughput drug screening in the PDX model and, based on these data, suggested a combinatorial treatment regimen targeting both EGFR and Src as the most efficient therapeutic option. In recognition of that drug responses are diverse in the pathologically identical tumors and even in the genetically homogeneous cancer cells [16, 17], profiling tumor transcriptome at the single-cell resolution can enable us to dissect inherent complex heterogeneous cell populations and to discern which cell would be different in drug responses affecting the ultimate outcomes that are covered under stochastic averaged signals in the bulk cell measurement. This strategy can be applied to other investigational models and eventually to patients with refractory cancer. Since higher engraftment rates are associated with more clinically aggressive tumors such as metastatic or treatment-refractory cases , surgically removed metastases or biopsy metastasis samples might extend the application of PDXs to advanced RCCs in terms of technical feasibility and unmet clinical needs.
Clonal genetic events in the metastases can be demonstrated for restricted subclones of the primary tumor, suggesting that only rare cells within the primary tumor have the ability to metastasize . The parallel progression model proposes independent evolution of early disseminates in distant sites, eventually leading to significant divergence between the primary and metastatic tumors . However, it remains largely unknown whether mRCCs also have initial ITH or whether the development of mRCCs is due to early dissemination or late diagnosis. Integrative analyses of PDXs derived from pRCC and distant metastasis to the lung from a single patient allowed us to infer which clonal and subclonal alterations contributed to tumor progression and spread to the metastatic lesion. The genetic differences between the pRCC and mRCC suggest highly complex non-linear branching clonal evolution as the basis for molecular and functional diversity, consistent with previous studies [7, 8]. Moreover, key pathways that are targeted by clinically available drugs showed distinct expression patterns between pRCC and mRCC. Most driver aberrations and actionable driver mutations that have therapeutic implications were located on the branches, suggesting that distinct intratumoral subclones had acquired different functional characteristics . The heterogeneity of metastatic cancer may underlie its poor responsiveness to therapy and explain why biomarkers of prognosis or therapy responsiveness measured exclusively from primary tumors give a restricted view of the biological properties of metastatic cancer .
Acquired resistance to therapy and disease progression can be due in part to intrinsic heterogeneity, including genetic, epigenetic, and biological properties of cancer cells that contribute to oncogenic activity. Thus, instead of single-agent therapy a rational approach that targets multiple subpopulations of tumor cells with a combination of non-cross-resistant drugs characterized by different mechanisms of action and non-overlapping profiles of toxicity will be necessary for long-lasting inhibition of highly heterogeneous mRCC [7, 8]. Knowledge of molecular alterations and features of tumors and the identification of mechanisms of tumor resistance provide the opportunity to test novel rationally designed drug combinations. Recent technological advances in single-cell sequencing have facilitated the paradigm shift in our understanding of the cancer ecosystem from the averaged signal of a complex tumor mass to the sum of distinctive signals in individual cells . Unique cellular behaviors reflecting oncogenic signatures have been extensively scrutinized by profiling transcriptomes from individual cells in various types of cancers [24, 47–49]. In the light of the precedent approaches for deconvolving heterogeneous cell populations and identifying specific cell types that are generally masked in bulk cell profiling, we could understand better the biologically relevant composition of cancer cells and their functional modulation within the tumor. Importantly, in this study such dissected intratumoral landscape enabled us to elicit the most effective drug combination of putting all the potentially targetable cancer cells together to be killed.
Using scRNA-seq, we could examine the heterogeneous drug target pathway activations at the single-cell level in a refractory mRCC patient. Distinct features of intratumoral expression variability across mRCC single cells that were masked in the bulk measurement prompted us to test the co-targeting strategy for the most vulnerable two signaling pathways with increased likelihood of complete response. Indeed, we observed significantly better treatment effects of the targeted combination therapy on mRCC-derived xenograft platforms in vitro and in vivo than monotherapies. Our findings described here will have a profound impact on translational research to overcome ITH-derived resistance and avoid ineffectual or unnecessary treatments. Although we could not analyze clinical response to our combination strategy because of rapid deterioration of the patient, we stress the utility and validity of single-cell transcriptome profiling in patients with refractory cancer for the design of personalized therapeutic strategies. In summary, the realization of the advantages in dissecting heterogeneous drug target pathway activations by scRNA-seq analysis will have significant clinical utility for the design of tailored combination therapy against highly heterogeneous tumors.
This study was carried out in accordance with the principles of the Declaration of Helsinki, and was approved by The Samsung Medical Center (Seoul, Korea) Institutional Review Board (IRB) (no. 2010-04-004). Participants in this study gave written informed consent for the research and publication of the results. Animal experiments were conducted in accordance with the Institute for Laboratory Animal Research Guide for the Care and Use of Laboratory Animals and the following protocols were approved by the IRB at the Samsung Medical Center (Seoul, Korea) (No. 20131217002). Animal care and handling was conducted in accordance with the National Institute of Health Guide for the Care and Use of Laboratory Animals (NIH publication no.80-23, revised 1978).
Sample selection and clinical characteristics of the patient
PDXs were established using surgically resected matched pRCC (pT1Nx; Fuhrman Grade 3) and paired lung metastasis from a 43-year-old man with ccRCC who experienced solitary lung metastasis 1 year after radical nephrectomy (IRB number 2010-04-004). He showed signs of rapid tumor dissemination to the bone, lung, pleura, and brain despite multiple salvage regimens, including pazopanib, everolimus, and high-dose interleukin (IL)-2, and finally died 16 months after complete metastatectomy as a result of rapid tumor progression. Serial 5-μm sections from each formalin-fixed paraffin-embedded block were processed for hematoxylin and eosin (H&E) staining and examined by a specialized pathologist.
Establishment of orthotopic PDXs
Athymic nude mice were obtained from Orient Bio (Seoul, Korea). Fresh tumor tissue was obtained from pRCC tissue of the patient by surgical excision under sterile conditions and matched fresh tumor tissue was taken from the mRCC. Each biopsied parental tumor mass was chopped into fragments, and was frozen or placed in formalin and embedded in paraffin for later analyses. A blood pellet was used for extraction of germline DNA. Fresh tumor tissue was stored on ice in Hank’s Balanced Salt Solution (Gibco, Grand Island, NY, USA) supplemented with penicillin/streptomycin (Gibco) for transport. For transplantation, 6- to 8-week-old NOD scid gamma mice were anesthetized with 100 mg/kg ketamine and 10 mg/kg xylazine. Primary tumor and paired lung metastasis samples minced into approximately 1-mm3 fragments in high-concentration Matrigel TM Basement Membrane Matrix (BD Biosciences, Franklin Lakes, NJ, USA) were directly implanted into the subrenal capsule (n = 4–5 for each tumor sample). For subrenal capsule implantation, xenograft tumor engraftment was defined as a palpable mass of >1 mm in diameter with pathologic confirmation. When the resulting tumors grew, each tumor (F1 generation) was resected as the primary tumor, divided, and passaged into five mice (F2 generation). The PDX tumors were harvested and divided into three samples for generation of second in vivo passage xenograft tumors, DNA/RNA extraction, and histopathologic examination. The origin of each xenograft was validated by short tandem repeat DNA fingerprinting. The process was repeated to produce subsequent generations via subcutaneous implantation in BALB/c nude mice to expand xenograft numbers.
We extracted genomic DNA from patient-derived tumor samples using the QIAamp DNA Minikit (Qiagen, Hilden, Germany) and from matched blood using the QIAamp DNA Blood kit (Qiagen). Purified DNA was sheared to an average size of 150 bp in a Covaris Adaptive Focused Acoustics™ (AFA) sonication device (S2, Covaris, Inc., Woburn, MA, USA) and indexed with unique barcode tags using PCR. Prepared libraries were assessed for quality and quantity using a Qubit 2.0 Fluorometer (Life Technologies, Carlsbad, CA, USA), 2100 Bioanalyzer (Agilent Inc., Palo Alto, CA, USA), and Mx3005P qPCR (Agilent Technologies, Inc). Exomes were targeted using the SureSelect XT Human All Exon V5 kit (Agilent). Samples were multiplexed and flow-cell clusters were created using the TruSeq Rapid Cluster kit and TruSeq Rapid SBS kit (Illumina, San Diego, CA, USA). Captured exomes were sequenced using the Illumina HiSeq 2500 platform, and paired-end 100-bp sequence data were generated. Sequencing reads that only mapped to the mouse genome reference (mm10) were filtered out.
Array comparative genomic hybridization
Purified DNA from patient-derived tumor samples was labeled with Cy5-dUTP following the Agilent Oligonucleotide Array-Based CGH for Genomic DNA Analysis protocol (Ver-7.3, Agilent). Cy5-labeled DNA was quantified together with reference DNA samples labeled with Cy3-dUTP to determine the DNA concentration using an ND-1000 Spectrophotometer (NanoDrop, Wilmington, DE, USA). Labeled test and reference samples were then hybridized to SurePrint G3 Human CGH 4 × 180 K Microarrays (Agilent) according to the manufacturer’s standard protocol. The dual-colored fluorescence signals were scanned using the Agilent Microarray Scanner and translated to log10 ratios using Feature Extraction software (Ver-18.104.22.168, Agilent).
Identification of single-nucleotide variants and copy number variants
Using the WES data, generated reads were aligned to the human genome reference hg19 using the Burrows-Wheeler Aligner  after removing duplicate reads, followed by implementation of the data-processing modules using the Genome Analysis Toolkit . Somatic SNVs were identified by Bayesian statistical analysis of bases and their qualities in the given tumor and paired normal BAM files at the genomic locus using the MuTect algorithm . Called variants were reviewed manually using the Integrative Genomic Viewer , and considered only within chromosomes 1–22 and X. Mutations were annotated using the SnpEff  package.
Based on the CGH data, extracted signals were normalized to log2 ratios using the limma package . To detect significant breakpoints across thousands of probe-derived signals, we applied the circular binary segmentation (CBS) algorithm using the DNAcopy package . After smoothing the data to detect outliers within autosomal chromosomes, aberrant segments were determined applying the significance level of 1.0E-04 to accept change-points based on a maximum t-statistic. We classified the segmented results into copy losses when the log2 ratios were lower than −0.25 and copy gains when these were greater than 0.25. Considering sample-specific tumor purity and ploidy, somatic copy number alterations (SCNA) were adjusted by implementing the ABSOLUTE algorithm . To compare the SCNA patterns across samples, segment values were averaged with 1-kb binning along the chromosomes.
Clustering of genomic clones
To determine the subclonal structure between primary RCC and paired lung metastasis, we adopted the PyClone algorithm  that computes the cellular prevalence of mutations and clusters these mutations based on a hierarchical Bayes statistical model. Mutational information, including somatic SNVs called from the deep exomes and absolute copy-number changes corresponding to SNV regions, was prepared for use with PyClone. The cellular prevalence for each mutation was estimated using a beta-binomial model by setting the number of Markov chain Monte Carlo (MCMC) iterations to 100,000, with a burn-in of 50,000. The number of clusters was inferred based on the average linkage hierarchical clustering in the post-burn-in trace by optimizing the maximization of posterior expected adjusted Rand index criterion.
Single-cell RNA sequencing and processing
scRNA-seq and data processing were carried out as previously described . Briefly, cells isolated from human mRCC, PDX-mRCC, and paired PDX-pRCC tumors were subjected to single-cell capture and cDNA amplification using the C1™ Single-Cell Auto Prep System (Fluidigm, South San Francisco, CA, USA) with the SMARTer kit (Clontech, Mountain View, CA, USA). RNA-seq reads for single cells and bulk samples were generated using the HiSeq 2500 in the 100-bp paired-end mode and reads that only mapped to the mouse genome reference (mm10) were subsequently removed. Filtered reads were aligned to the human genome reference hg19 with the sample specific-splice junction using the 2-pass mode of STAR (Ver-2.4.0d) . Transcripts per million (TPM) was quantified using RSEM (Ver-1.2.18) . To filter out poor quality cells, we applied the critea of >1 M reads per cell, >60 % uniquely mapped rate, >35 % exonic region coverage rate, and >5000 detected genes. We considered TPM values greater than 1 to be reliable and only focused on genes that were detected in more than 10 % of a group between pRCC and paired lung metastasis. To identify differential and common expression signatures between pRCC and paired lung metastasis single cells compared to normal signals, we used expression profiles of normal kidney cortex from the GTEx portal (http://www.gtexportal.org/home/; transcript read counts V3) by converting to TPM values. Finally, we normalized sample-to-sample variation by applying a mean centroid.
Estimation of activity for expression signatures and drug sensitivities
To understand the relative activation status for a pathway or signature across samples, we implemented the GSVA algorithm  in RNA-seq data. GSVA scores were estimated for given gene sets from the following sources: stromal signature, extracted from the ESTIMATE package ; ccRCC signature, taken from Jones et al. ; EMT-induced signature, taken from Taube et al. ; metastatic signature, taken from Jones et al. ; prognostic signature, taken from The Cancer Genome Atlas Research Network for ccRCC ; EGFR signaling, (Reactome, signaling by constitutively active EGFR); Src signaling, extracted from Gatza M.L. et al. ; mTOR signaling (Reactome, mTOR signaling); VEGFR signaling (PID, signaling events mediated by VEGFR1 and VEGFR2); RAF signaling, (Reactome, RAF activation); MEK signaling, (Reactome, MEK activation); c-Met signaling, (PID, signaling events mediated by c-Met); SCF-KIT signaling, (Reactome, signaling by SCF-KIT); PI3K/AKT signaling (Reactome, PI3K/AKT signaling in cancer); FGFR signaling (Reactome, signaling by FGFR); PDGFR signaling (PID, PDGF receptor signaling network). To evaluate whether an estimated gene set signature is significantly activated, we transformed the observed GSVA scores to binary scores. Gene sets with the same size as each original panel of genes were randomly generated with permutation (×1000) and computed for their GSVA scores. Assignment of the original GSVA scores as “activated” was determined with the cutoff values of averaged scores in the randomly selected gene sets.
In addition to estimation of activity status for signaling pathways, relevant targeting drug sensitivities were also predicted from expression profiles. Following a previous approach , cancer cell line expression data with measured drug response data from the Cancer Genome Project (CGP)  were used as a training set. After adjusting two independent sets of our expression data with the training set of solid tumor cell lines using ComBat , a ridge regression model was fitted for the training set with the given drugs that were simultaneously identified in our study. The lowest varying 20 % of genes were filtered out to focus on biological variability over technical variability. To evaluate the prediction sensitivity, leave-one-out-cross-validation (LOOCV) was applied with the total dataset. From this computation for a total of 10 drugs (BIBW2992 [Afatinib], Dasatinib, Gefitinib, Erlotinib, Temsirolimus, Pazopanib, Sunitinib, Sorafenib, AZD6244 [Selumetinib], and PF-02341066 [Crizotinib]), our data were tested to predict drug sensitivity and the distribution of estimates was transformed to Z-scores by dividing the averaged drug sensitivity estimates by the standard deviation for the difference. To determine the prediction accuracy of drug sensitivity estimates with measured drug sensitivity from high-throughput drug screening, nanomolar scaled IC50 values were also transformed to Z-scores.
Functional network analysis
To identify discrete biological networks respectively enriched in pRCC and mRCC cells, differentially expressed genes (DEGs) were first defined with a criteria of the Benjamini–Hochberg corrected FDR <0.01 and fold changes of at least twofold. DEGs were separately applied to Gene Ontology (GO) category analysis using the ClueGO  plug-in within the Cytoscape framework . The statistical significance for the over-represented pathways in the GO Biological Process category was estimated using Benjamini–Hochberg correction for multiple testing-controlled P values. Significantly enriched terms were functionally grouped based on kappa scores >0.3.
Primary in vitro short-term culture
Xenograft tumor specimens were dissociated into single cells according to previously published protocols . Dissociated PDX cells were cultured in neurobasal media-A supplemented with N2 (×1/2, Life Technologies, Carlsbad, CA, USA), B27 (×1/2, Gibco), basic fibroblast growth factor (20 ng/mL; R&D Systems, Minneapolis, MN, USA), epidermal growth factor (EGF, 20 ng/mL; R&D Systems), neuregulin 1 (10 ng/mL; R&D Systems), and insulin-like growth factor 1 (100 ng/mL; R&D Systems) and containing 10 % conditioned medium (CM) from human mesenchymal stem cells (MSCs). To generate the CM, MSCs were grown to 70 % confluency in plates with 10 % FBS L-DMEM and allowed to adhere overnight at 37 °C and 5 % CO2. The following day, the medium was replaced with serum-free culture medium and the cells were cultured for 2 days. The used medium was collected as MSC-CM, centrifuged to remove cell debris, and passed through a 0.45-μm filter. CM aliquots were frozen at −80 °C until use.
In vitro drug sensitivity assay
Primary RCC PDX cells cultured under serum-free sphere culture conditions were seeded in 384-well plates at 500 cells per well. Two hours after plating, the cells were treated with a drug library in threefold and 10-point serial dilution series (n = 3 for each condition). After incubation at 37 °C in a 5 % CO2 humidified incubator for 6 days, cell viability was analyzed using an adenosine triphosphate monitoring system based on firefly luciferase (ATPliteTM 1step, PerkinElmer, Waltham, MA, USA). The drug library was composed of 20 targeted agents that were included in the clinical guidelines or in current clinical trials (gefitinib, erlotinib, lapatinib, afatinib, tivantinib, foretinib, crizotinib, selumetinib, temsirolimus, everolimus, cabozantinib, vandetanib, sunitinib, sorafenib, dovitinib, vemurafenib, BKM 120, pazopanib, nintedanib, and DAPT; all purchased from Selleckchem, Houston, TX, USA). The drugs were stored and diluted according to the manufacturer’s instructions. Test concentrations for each drug were empirically derived to produce a clinically relevant spectrum of drug activity. Half-maximal (50 %) inhibitory concentration values (IC50) were calculated as an average of triplicate experiments using the S+ Chip Analyzer (Samsung Electro-Mechanics Company, Ltd., Gyeonggi, Korea) .
For signal transduction assays under treatment with the targeted agents, primary cultured PDX cells were maintained overnight in serum-free sphere culture conditions without growth factors, incubated for 1 h with each inhibitor, and pulsed with original culture medium supplemented for 15 min. For immunoblotting, cells were lysed in RIPA lysis buffer supplemented with 1× phosphatase inhibitors (PhosStop; Roche Diagnostics, Basel, Switzerland) and a 1× protease inhibitor cocktail (Complete Mini; Roche Diagnostics). After centrifugation at 10,000 × g for 5 min, the supernatant was harvested and protein concentration was determined using a bicinchoninic acid protein assay kit (Thermo Scientific, Waltham, MA, USA). Equal amounts of protein were subjected to SDS-PAGE and transferred to polyvinylidene difluoride membranes (Whatman plc, Little Chalfont, UK). Membranes were blocked in 5 % skim milk or bovine serum albumin for 1 h at room temperature, incubated with the indicated primary antibodies overnight, and then with the appropriate secondary antibodies. Antibodies against p-EGFR (Tyr1068), p-Src (Tyr527), p-ERK, p-AKT (Ser473) (all purchased from Cell Signaling Technology, Danvers, MA, USA), and GAPDH (Santa Cruz Biotechnology, Santa Cruz, CA, USA) were used.
Preparation of microfluidic drug screening device
The microfluidic drug screening device was made of polydimethylsiloxane (PDMS, Sylgard 184; Dow Corning, Corning, NY, USA) using a conventional softlithography process with 200 micron high SU-8 patterned silicon wafer (MicroChem, Westborough, MA, USA). The fabricated device was sterilized and bonded onto a cover glass to enclose the microchannels by oxygen plasma treatment (Femto Science, Somerset, NJ, USA). The device was then heated in an oven at 80 °C for 24 h to allow the surface to recover its hydrophobicity. The restored hydrophobicity of the microfluidic channel surface helps the injected extracellular matrix (ECM) form a stable interface with the side channels. In repeated measures analysis for drug sensitivity, the identical experimental condition was applied to the second measurement in PDX-mRCC cells.
Three-dimensional cell culture and drug treatment
Collagen type 1 (3 mg/mL, rat tail; Corning) was used as an ECM scaffold to embed the cells. The collagen solution was prepared at pH 7 and a concentration of 2 mg/mL. This solution was diluted in a mixture of 10× phosphate-buffered saline (PBS; Gibco) and sterilized deionized water. The pH was adjusted using 0.5 N NaOH. Dissociated cells were suspended in the collagen solution at a density of 0.5 × 106 cells/mL. The suspension was injected into a center channel of the device and allowed to gel by incubation at 37 °C and 5 % CO2 for 30 min. Details of the device preparation and gel filling procedure have been described previously . To avoid cell attachment to the microfluidic channel surface, the device was inverted every 5 min. After gel formation, the side channels were filled with medium containing each drug candidate. The medium in the channel was refreshed every 24 h.
Cell viability was quantified at 4 and 7 days of culture using the Live/Dead Viability Assay Kit (Molecular Probes, Invitrogen, Carlsbad, CA, USA) containing calcein AM and ethidium homodimer to identify live (green) and dead (red) cells, respectively. Cells in the microfluidic device were incubated at 37 °C with 5 % CO2 for 30 min, and then the staining solution was replaced with PBS. The cells were counted using ImageJ software (Image Processing and Analysis in Java, NIH, Bethesda, MD, USA). Cell viability was calculated as the number of live cells divided by the total cell number. Normalized cell viability was determined by dividing cell viability by the viability of cells cultured in pure medium condition.
In vivo drug efficacy
BALB/c nude mice (female, 6–8 weeks old; Orient Bio) were used for drug efficacy studies. Animal experiments were conducted in accordance with the Institute for Laboratory Animal Research Guide for the Care and Use of Laboratory Animals and the following protocols were approved by the IRB at Samsung Medical Center (Seoul, Korea). Briefly, 2 × 105 dissociated PDX-mRCC cells mixed 1:1 with Matrigel (BD Biosciences) were inoculated subcutaneously into the right flank of each mouse. Tumor diameters were measured with calipers twice a week and tumor volume in mm3 was calculated using the following formula: tumor volume = (l × w2)/2, where l is the longest diameter of the tumor and w is the shortest diameter of the tumor. Mice bearing established tumors (100–150 mm3) were randomly allocated to four groups (5 in each group): vehicle, afatinib at 20 mg/kg (p.o.), dasatinib 15 mg/kg (i.p.), or afatinib/dasatinib combination on a daily dosing regimen for up to 15 days. Mice in the four groups exhibited similar average tumor volumes and body weight. Throughout the study, the mice were weighed and tumor burden was monitored every 3 days. Mean tumor volumes were calculated, and growth curves were established as a function of time. The error bars indicate the value of the standard error of the mean (SEM). Tumors from each group were collected at the end of the experiment for further analysis.
Immunohistochemical and TUNEL staining
Tumors were embedded in paraffin, sectioned at 5 μm, and stained with H&E. Paraffin-embedded tissue sections were deparaffinized and rehydrated. For immunochemical staining, heat-induced epitope retrieval was performed using a target retrieval solution (Dako, Glostrup, Denmark) for 5 min in a microwave. Slides were treated with 3 % hydrogen peroxide for 10 min to inactivate endogenous peroxidase and then blocked for 20 min at room temperature (RT) in blocking solution (5 % normal horse serum, 1 % normal goat serum, 0.1 % Triton-X 100 in 1× PBS). After blocking, the slides were incubated with primary antibodies, including mouse monoclonal antibody against p-ERK (Cell Signaling Technology) and human Ki-67 (BD Pharmingen). After washing, the slides were incubated with secondary antibodies for 1 h at RT, and counterstained with hematoxylin (Sigma-Aldrich, St. Louis, MO, USA). Apoptotic cells were identified by histologic analysis of DNA fragmentation in paraffin sections of the xenograft tumors. We performed terminal deoxynucleotide transferase-mediated dUTP nick end labeling (TUNEL) staining on the tumor sections using the DeadEnd™ Colorimetric TUNEL System (Promega, Madison, WI, USA). The proliferative or apoptotic index was calculated as a ratio of Ki-67-positive or TUNEL-positive cell number to total cell number in high-power (×400) fields.
All values are expressed as the mean ± SEM. Comparisons between two groups were analyzed by Student’s t-test. One-way analysis of variance was applied for comparisons between more than two groups and to determine the statistical significance for the fitting model in the linear regression of two components. All P values were two-sided, and P <0.05 was considered statistically significant. For discovery of differentially expressed genes (DEGs) between pRCC and mRCC cells, we applied the Benjamini–Hochberg correction for multiple-testing. Statistically significant DEGs were regarded using the cutoff of FDR <0.01 and fold-change ≥2. Multiple regression analysis was applied to evaluate the transcriptomic heterogeneity across single cells, and to estimate the explanatory power of randomly selected single cells (with the given number by permutation × 1000) attributed to the expression profile of the bulk measurement by calculating the adjusted R-square. All data analyses were performed using SPSS statistical software, version 19.0 (SPSS, Inc., Chicago, IL, USA).
RNA sequence and aCGH data have been deposited in the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) under accession number GSE73122. Whole-exome sequence data can be accessed at the NCBI Sequence Read Archive (SRA) with accession number SRP063388.
We thank research fellows in the Samsung Genome Institute for experimental assistance and valuable comments on our work. This research was supported by the Samsung Medical Center and Health Technology R&D Project, Ministry of Health & Welfare (HI13C2096000013).
Open AccessThis 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.
- Rini BI, Atkins MB. Resistance to targeted therapy in renal-cell carcinoma. Lancet Oncol. 2009;10:992–1000.View ArticlePubMedGoogle Scholar
- Linehan WM, Srinivasan R, Schmidt LS. The genetic basis of kidney cancer: a metabolic disease. Nat Rev Urol. 2010;7:277–85.View ArticlePubMedPubMed CentralGoogle Scholar
- Ravaud A, Hawkins R, Gardner JP, von der Maase H, Zantl N, Harper P, Rolland F, Audhuy B, Machiels JP, Petavy F, et al. Lapatinib versus hormone therapy in patients with advanced renal cell carcinoma: a randomized phase III clinical trial. J Clin Oncol. 2008;26:2285–91.Google Scholar
- Rini BI, Escudier B, Tomczak P, Kaprin A, Szczylik C, Hutson TE, Michaelson MD, Gorbunova VA, Gore ME, Rusakov IG, et al. Comparative effectiveness of axitinib versus sorafenib in advanced renal cell carcinoma (AXIS): a randomised phase 3 trial. Lancet. 2011;378:1931–9.Google Scholar
- Hudes G, Carducci M, Tomczak P, Dutcher J, Figlin R, Kapoor A, Staroslawska E, Sosman J, McDermott D, Bodrogi I, et al. Temsirolimus, interferon alfa, or both for advanced renal-cell carcinoma. N Engl J Med. 2007;356:2271–81.Google Scholar
- Gupta K, Miller JD, Li JZ, Russell MW, Charbonneau C. Epidemiologic and socioeconomic burden of metastatic renal cell carcinoma (mRCC): a literature review. Cancer Treat Rev. 2008;34:193–205.View ArticlePubMedGoogle Scholar
- Gerlinger M, Rowan AJ, Horswell S, Larkin J, Endesfelder D, Gronroos E, Martinez P, Matthews N, Stewart A, Tarpey P, et al. Intratumor heterogeneity and branched evolution revealed by multiregion sequencing. N Engl J Med. 2012;366:883–92.Google Scholar
- Gerlinger M, Horswell S, Larkin J, Rowan AJ, Salm MP, Varela I, Fisher R, McGranahan N, Matthews N, Santos CR, et al. Genomic architecture and evolution of clear cell renal cell carcinomas defined by multiregion sequencing. Nat Genet. 2014;46:225–33.Google Scholar
- Kosari F, Parker AS, Kube DM, Lohse CM, Leibovich BC, Blute ML, Cheville JC, Vasmatzis G. Clear cell renal cell carcinoma: gene expression analyses identify a potential signature for tumor aggressiveness. Clin Cancer Res. 2005;11:5128–39.Google Scholar
- Ozcan A, de la Roza G, Ro JY, Shen SS, Truong LD. PAX2 and PAX8 expression in primary and metastatic renal tumors: a comprehensive comparison. Arch Pathol Lab Med. 2012;136:1541–51.View ArticlePubMedGoogle Scholar
- Wyler L, Napoli CU, Ingold B, Sulser T, Heikenwalder M, Schraml P, Moch H. Brain metastasis in renal cancer patients: metastatic pattern, tumour-associated macrophages and chemokine/chemoreceptor expression. Br J Cancer. 2014;110:686–94.Google Scholar
- Burrell RA, McGranahan N, Bartek J, Swanton C. The causes and consequences of genetic heterogeneity in cancer evolution. Nature. 2013;501:338–45.View ArticlePubMedGoogle Scholar
- Higashiyama M, Okami J, Maeda J, Tokunaga T, Fujiwara A, Kodama K, Imamura F, Kobayashi H. Differences in chemosensitivity between primary and paired metastatic lung cancer tissues: In vitro analysis based on the collagen gel droplet embedded culture drug test (CD-DST). J Thorac Dis. 2012;4:40–7.Google Scholar
- Xu X, Hou Y, Yin X, Bao L, Tang A, Song L, Li F, Tsang S, Wu K, Wu H, et al. Single-cell exome sequencing reveals single-nucleotide mutation characteristics of a kidney tumor. Cell. 2012;148:886–95.Google Scholar
- Dancey JE, Bedard PL, Onetto N, Hudson TJ. The genetic basis for cancer treatment decisions. Cell. 2012;148:409–20.View ArticlePubMedGoogle Scholar
- Cohen AA, Geva-Zatorsky N, Eden E, Frenkel-Morgenstern M, Issaeva I, Sigal A, Milo R, Cohen-Saidon C, Liron Y, Kam Z, et al. Dynamic proteomics of individual cancer cells in response to a drug. Science. 2008;322:1511–6.Google Scholar
- Spencer SL, Gaudet S, Albeck JG, Burke JM, Sorger PK. Non-genetic origins of cell-to-cell variability in TRAIL-induced apoptosis. Nature. 2009;459:428–32.View ArticlePubMedPubMed CentralGoogle Scholar
- Costello JC, Heiser LM, Georgii E, Gonen M, Menden MP, Wang NJ, Bansal M, Ammadud-din M, Hintsanen P, Khan SA, et al et al. A community effort to assess and improve drug sensitivity prediction algorithms. Nat Biotechnol. 2014;32:1202–12.Google Scholar
- Klijn C, Durinck S, Stawiski EW, Haverty PM, Jiang Z, Liu H, Degenhardt J, Mayba O, Gnad F, Liu J, et al. A comprehensive transcriptional portrait of human cancer cell lines. Nat Biotechnol. 2015;33:306–12.Google Scholar
- Julien S, Merino-Trigo A, Lacroix L, Pocard M, Goere D, Mariani P, Landron S, Bigot L, Nemati F, Dartigues P, et al. Characterization of a large panel of patient-derived tumor xenografts representing the clinical heterogeneity of human colorectal cancer. Clin Cancer Res. 2012;18:5314–28.Google Scholar
- Ding L, Ellis MJ, Li S, Larson DE, Chen K, Wallis JW, Harris CC, McLellan MD, Fulton RS, Fulton LL, et al. Genome remodelling in a basal-like breast cancer metastasis and xenograft. Nature. 2010;464:999–1005.Google Scholar
- DeRose YS, Wang G, Lin YC, Bernard PS, Buys SS, Ebbert MT, Factor R, Matsen C, Milash BA, Nelson E, et al. Tumor grafts derived from women with breast cancer authentically reflect tumor pathology, growth, metastasis and disease outcomes. Nat Med. 2011;17:1514–20.Google Scholar
- Kreso A, O’Brien CA, van Galen P, Gan OI, Notta F, Brown AM, Ng K, Ma J, Wienholds E, Dunant C, et al. Variable clonal repopulation dynamics influence chemotherapy response in colorectal cancer. Science. 2013;339:543–8.Google Scholar
- Kim KT, Lee HW, Lee HO, Kim SC, Seo YJ, Chung W, Eum HH, Nam DH, Kim J, Joo KM, Park WY. Single-cell mRNA sequencing identifies subclonal heterogeneity in anti-cancer drug responses of lung adenocarcinoma cells. Genome Biol. 2015;16:127.Google Scholar
- The Cancer Genome Atlas Research Network. Comprehensive molecular characterization of clear cell renal cell carcinoma. Nature. 2013;499:43–9.View ArticlePubMed CentralGoogle Scholar
- Erdogan F, Demirel A, Polat O. Prognostic significance of morphologic parameters in renal cell carcinoma. Int J Clin Pract. 2004;58:333–6.View ArticlePubMedGoogle Scholar
- Tsai JH, Yang J. Epithelial-mesenchymal plasticity in carcinoma metastasis. Genes Dev. 2013;27:2192–206.View ArticlePubMedPubMed CentralGoogle Scholar
- Mikami S, Oya M, Mizuno R, Kosaka T, Katsube K, Okada Y. Invasion and metastasis of renal cell carcinoma. Med Mol Morphol. 2014;47:63–7.View ArticlePubMedGoogle Scholar
- Heng DY, Xie W, Regan MM, Harshman LC, Bjarnason GA, Vaishampayan UN, Mackenzie M, Wood L, Donskov F, Tan MH, et al. External validation and comparison with other models of the International Metastatic Renal-Cell Carcinoma Database Consortium prognostic model: a population-based study. Lancet Oncol. 2013;14:141–8.Google Scholar
- Geeleher P, Cox NJ, Huang RS. Clinical drug response can be predicted using baseline gene expression levels and in vitro drug sensitivity in cell lines. Genome Biol. 2014;15:R47.View ArticlePubMedPubMed CentralGoogle Scholar
- Garnett MJ, Edelman EJ, Heidorn SJ, Greenman CD, Dastur A, Lau KW, Greninger P, Thompson IR, Luo X, Soares J, et al. Systematic identification of genomic markers of drug sensitivity in cancer cells. Nature. 2012;483:570–5.Google Scholar
- Smith K, Gunaratnam L, Morley M, Franovic A, Mekhail K, Lee S. Silencing of epidermal growth factor receptor suppresses hypoxia-inducible factor-2-driven VHL−/− renal cancer. Cancer Res. 2005;65:5221–30.View ArticlePubMedGoogle Scholar
- Moch H, Sauter G, Buchholz N, Gasser TC, Bubendorf L, Waldman FM, Mihatsch MJ. Epidermal growth factor receptor expression is associated with rapid tumor cell proliferation in renal cell carcinoma. Hum Pathol. 1997;28:1255–9.Google Scholar
- Stumm G, Eberwein S, Rostock-Wolf S, Stein H, Pomer S, Schlegel J, Waldherr R. Concomitant overexpression of the EGFR and erbB-2 genes in renal cell carcinoma (RCC) is correlated with dedifferentiation and metastasis. Int J Cancer. 1996;69:17–22.Google Scholar
- Yonezawa Y, Nagashima Y, Sato H, Virgona N, Fukumoto K, Shirai S, Hagiwara H, Seki T, Ariga T, Senba H, et al. Contribution of the Src family of kinases to the appearance of malignant phenotypes in renal cancer cells. Mol Carcinog. 2005;43:188–97.Google Scholar
- Zhao Y, Li X, Sun X, Zhang Y, Ren H. EMT phenotype is induced by increased Src kinase activity via Src-mediated caspase-8 phosphorylation. Cell Physiol Biochem. 2012;29:341–52.View ArticlePubMedGoogle Scholar
- Keating GM. Afatinib: a review of its use in the treatment of advanced non-small cell lung cancer. Drugs. 2014;74:207–21.View ArticlePubMedGoogle Scholar
- Chang Q, Jorgensen C, Pawson T, Hedley DW. Effects of dasatinib on EphA2 receptor tyrosine kinase activity and downstream signalling in pancreatic cancer. Br J Cancer. 2008;99:1074–82.View ArticlePubMedPubMed CentralGoogle Scholar
- Montero JC, Seoane S, Ocana A, Pandiella A. Inhibition of SRC family kinases and receptor tyrosine kinases by dasatinib: possible combinations in solid tumors. Clin Cancer Res. 2011;17:5546–52.View ArticlePubMedGoogle Scholar
- Friedrich J, Seidel C, Ebner R, Kunz-Schughart LA. Spheroid-based drug screen: considerations and practical approach. Nat Protoc. 2009;4:309–24.View ArticlePubMedGoogle Scholar
- Shin Y, Han S, Jeon JS, Yamamoto K, Zervantonakis IK, Sudo R, Kamm RD, Chung S. Microfluidic assay for simultaneous culture of multiple cell types on surfaces or within hydrogels. Nat Protoc. 2012;7:1247–59.Google Scholar
- Navin NE. Cancer genomics: one cell at a time. Genome Biol. 2014;15:452.View ArticlePubMedPubMed CentralGoogle Scholar
- Siolas D, Hannon GJ. Patient-derived tumor xenografts: transforming clinical samples into mouse models. Cancer Res. 2013;73:5315–9.View ArticlePubMedPubMed CentralGoogle Scholar
- Wu X, Northcott PA, Dubuc A, Dupuy AJ, Shih DJ, Witt H, Croul S, Bouffet E, Fults DW, Eberhart CG, et al. Clonal selection drives genetic divergence of metastatic medulloblastoma. Nature. 2012;482:529–33.Google Scholar
- Klein CA. Parallel progression of primary tumours and metastases. Nat Rev Cancer. 2009;9:302–12.View ArticlePubMedGoogle Scholar
- Yokota J, Kohno T. Molecular footprints of human lung cancer progression. Cancer Sci. 2004;95:197–204.View ArticlePubMedGoogle Scholar
- Patel AP, Tirosh I, Trombetta JJ, Shalek AK, Gillespie SM, Wakimoto H, Cahill DP, Nahed BV, Curry WT, Martuza RL, et al. Single-cell RNA-seq highlights intratumoral heterogeneity in primary glioblastoma. Science. 2014;344:1396–401.Google Scholar
- Lee MC, Lopez-Diaz FJ, Khan SY, Tariq MA, Dayn Y, Vaske CJ, Radenbaugh AJ, Kim HJ, Emerson BM, Pourmand N. Single-cell analyses of transcriptional heterogeneity during drug tolerance transition in cancer cells by RNA sequencing. Proc Natl Acad Sci U S A. 2014;111:E4726–35.Google Scholar
- Yu M, Ting DT, Stott SL, Wittner BS, Ozsolak F, Paul S, Ciciliano JC, Smas ME, Winokur D, Gilman AJ, et al. RNA sequencing of pancreatic circulating tumour cells implicates WNT signalling in metastasis. Nature. 2012;487:510–3.Google Scholar
- Li H, Durbin R. Fast and accurate long-read alignment with Burrows-Wheeler transform. Bioinformatics. 2010;26:589–95.View ArticlePubMedPubMed CentralGoogle Scholar
- McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, Garimella K, Altshuler D, Gabriel S, Daly M, DePristo MA. The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 2010;20:1297–303.Google Scholar
- Cibulskis K, Lawrence MS, Carter SL, Sivachenko A, Jaffe D, Sougnez C, Gabriel S, Meyerson M, Lander ES, Getz G. Sensitive detection of somatic point mutations in impure and heterogeneous cancer samples. Nat Biotechnol. 2013;31:213–9.Google Scholar
- Robinson JT, Thorvaldsdottir H, Winckler W, Guttman M, Lander ES, Getz G, Mesirov JP. Integrative genomics viewer. Nat Biotechnol. 2011;29:24–6.Google Scholar
- Cingolani P, Platts A, le Wang L, Coon M, Nguyen T, Wang L, Land SJ, Lu X, Ruden DM. A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3. Fly (Austin). 2012;6:80–92.Google Scholar
- Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43:e47.Google Scholar
- Venkatraman ES, Olshen AB. A faster circular binary segmentation algorithm for the analysis of array CGH data. Bioinformatics. 2007;23:657–63.View ArticlePubMedGoogle Scholar
- Carter SL, Cibulskis K, Helman E, McKenna A, Shen H, Zack T, Laird PW, Onofrio RC, Winckler W, Weir BA, et al. Absolute quantification of somatic DNA alterations in human cancer. Nat Biotechnol. 2012;30:413–21.Google Scholar
- Roth A, Khattra J, Yap D, Wan A, Laks E, Biele J, Ha G, Aparicio S, Bouchard-Cote A, Shah SP. PyClone: statistical inference of clonal population structure in cancer. Nat Methods. 2014;11:396–8.Google Scholar
- Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, Batut P, Chaisson M, Gingeras TR. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29:15–21.Google Scholar
- Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12:323.View ArticlePubMedPubMed CentralGoogle Scholar
- Hanzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013;14:7.View ArticlePubMedPubMed CentralGoogle Scholar
- Yoshihara K, Shahmoradgoli M, Martinez E, Vegesna R, Kim H, Torres-Garcia W, Trevino V, Shen H, Laird PW, Levine DA, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013;4:2612.Google Scholar
- Jones J, Otu H, Spentzos D, Kolia S, Inan M, Beecken WD, Fellbaum C, Gu X, Joseph M, Pantuck AJ, et al. Gene signatures of progression and metastasis in renal cell cancer. Clin Cancer Res. 2005;11:5730–9.Google Scholar
- Taube JH, Herschkowitz JI, Komurov K, Zhou AY, Gupta S, Yang J, Hartwell K, Onder TT, Gupta PB, Evans KW, et al. Core epithelial-to-mesenchymal transition interactome gene-expression signature is associated with claudin-low and metaplastic breast cancer subtypes. Proc Natl Acad Sci U S A. 2010;107:15449–54.Google Scholar
- Gatza ML, Lucas JE, Barry WT, Kim JW, Wang Q, Crawford MD, Datto MB, Kelley M, Mathey-Prevot B, Potti A, Nevins JR. A pathway-based classification of human breast cancer. Proc Natl Acad Sci U S A. 2010;107:6994–9.Google Scholar
- Johnson WE, Li C, Rabinovic A. Adjusting batch effects in microarray expression data using empirical Bayes methods. Biostatistics. 2007;8:118–27.View ArticlePubMedGoogle Scholar
- Bindea G, Mlecnik B, Hackl H, Charoentong P, Tosolini M, Kirilovsky A, Fridman WH, Pages F, Trajanoski Z, Galon J. ClueGO: a Cytoscape plug-in to decipher functionally grouped gene ontology and pathway annotation networks. Bioinformatics. 2009;25:1091–3.Google Scholar
- Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, Amin N, Schwikowski B, Ideker T. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13:2498–504.Google Scholar
- Lee HW, Lee JI, Lee SJ, Cho HJ, Song HJ, da Jeong E, Seo YJ, Shin S, Joung JG, Kwon YJ, et al. Patient-derived xenografts from non-small cell lung cancer brain metastases are valuable translational platforms for the development of personalized targeted therapy. Clin Cancer Res. 2015;21:1172–82.Google Scholar
- Lee DW, Choi YS, Seo YJ, Lee MY, Jeon SY, Ku B, Kim S, Yi SH, Nam DH. High-throughput screening (HTS) of anticancer drug efficacy on a micropillar/microwell chip platform. Anal Chem. 2014;86:535–42.Google Scholar