Precision medicine for hepatocelluar carcinoma using molecular pattern diagnostics: results from a preclinical pilot study

The aim of this study was to design a road map for personalizing cancer therapy in hepatocellular carcinoma (HCC) by using molecular pattern diagnostics. As an exploratory study, we investigated molecular patterns of tissues of two tumors from individual HCC patients, which in previous experiments had shown contrasting reactions to the phase 2 transforming growth factor beta receptor 1 inhibitor galunisertib. Cancer-driving molecular patterns encompass – inter alias – altered transcription profiles and somatic mutations in coding regions differentiating tumors from their respective peritumoral tissues and from each other. Massive analysis of cDNA ends and all-exome sequencing demonstrate a highly divergent transcriptional and mutational landscape, respectively, for the two tumors, that offers potential explanations for the tumors contrasting responses to galunisertib. Molecular pattern diagnostics (MPDs) suggest alternative, individual-tumor-specific therapies, which in both cases deviate from the standard sorafenib treatment and from each other. Suggested personalized therapies use kinase inhibitors and immune-focused drugs as well as low-toxicity natural compounds identified using an advanced bioinformatics routine included in the MPD protocol. The MPD pipeline we describe here for the prediction of suitable drugs for treatment of two contrasting HCCs may serve as a blueprint for the design of therapies for various types of cancer.

Hepatocellular carcinoma (HCC) is one of the most lethal cancers worldwide. Nearly 745 000 people died from it only in 2012. 1 Patients' 5-year overall survival (OS) rate of o20% indicates the urgent need for alternative therapies to improve the outcome for these patients. 2 HCC develops along different clinical histories including chronic hepatitis, cirrhosis and alcoholism. 3 All these factors contribute to unceasing inflammation and regeneration of hepatocytes, making it challenging to achieve diagnosis and prognosis of HCC at earlier stages.
Currently, the multikinase inhibitor sorafenib is the only effective, approved systemic therapy for advanced HCC that are not suitable for other curative treatment, 4 but the occurrence of side effects has markedly reduced the impact of the drug in daily life clinical practice. 5,6 Given the limited efficiency of the standard treatment, as well as the occurrence of drug resistance, 7 we addressed the question whether the newly arising concept of precision oncology could enable us to design novel therapeutic strategies that take into account the genetic diversity of these patients' tumors.
Other drugs such as the transforming growth factor beta receptor 1 (TGFBR1) blocker galunisertib are undergoing clinical trials for the treatment of HCC. 8 A large body of evidence indicates that TGFB1 is an important key to tumor progression, as it promotes the epithelial-to-mesenchymal transition (EMT) and activates the WNT pathway, a hallmark of HCC. 9 Immune therapy is currently being considered for the treatment of HCC and a comprehensive meta-analysis of recent studies encompassing more than 1800 patients indicates that patients undergoing specific immunotherapy benefit from a significantly higher overall and recurrence-free survival than those in control groups. 10 TGFB1 plays an important role in the regulation of immune responses via cancer-associated fibroblasts (CAFs) that express the growth factor in a self-sustaining autocrine cycle. CAFs sustain oncogenic features of cancer cells including suppression of the functions of various immune cells, particularly effector T cells and natural killer (NK) cells. TGFB1 also regulates T-regulatory cells (Treg) maturation and thereby suppresses immune responses. 11 As in other cancers there is ample evidence that also in HCC, control of the immune system by the neoplastic complex contributes significantly to the survival of cancer cells. It has been shown, for example, that the presence of a certain dysfunctional subset of tumor-infiltrating NK cells is associated with tumor progression and is an independent indicator of poor outcome in HCC patients. 12 Recent work indicates that patients suffering from refractory cancers that were treated by genomics-guided precision medicine did indeed have a significantly better progress-free survival (PFS) ratios and longer median PFS compared with patients who did not receive personalized therapy. 13 Precision oncology has been facilitated by the advent of next-generation sequencing, which enables particular molecular genetic profiles to be identified in the individual patient who may be targeted by precise, personalized therapy. Identified targets are then used to search databases for drugs that address these aberrantly expressed molecules and pathways using the bioinformatics pipeline. This concept benefits from the fact that drugs have been developed and are applied for many molecular targets, across a plethora of different diseases. The knowledge of the individual architecture of a patient's cancer may now enable these drugs to be directed against these specific oncogenic features, in a form of one-person trial. Drug repurposing, retasking or reprofiling has already been demonstrated as a promising strategy for cancer therapy, 14 which could be justified if the patient's tumor reveals molecular patterns indicative of an altered mitochondrial function, like the Warburg effect. Thus, given that suitable targets can be identified for the individual cancer, readdressable drugs are often at hand to attack it.
Recently, we have characterized HCC tissues by their differential response to galunisertib 15 using NGS-based massive analysis of cDNA ends (MACE), 16,17 highthroughput transcription profiling for the investigation of aberrant signaling and metabolic pathways and all-exomesequencing for the identification of somatic mutations accompanying neoplasia. 18 Based on this comprehensive set of data, molecular pattern diagnostics (MPDs) first identify biological pathways most enriched in overexpressed genes in the tumor. MPDs then help to select the potentially best suited drugs according to the number of targeted enriched pathways, resulting in a ranked list of drugs suited for approved use or repurposing in the frame of a monotherapy or for combination therapies with other immune-related drugs or natural compounds.
The aim of this study was to design a road map for precision medicine of HCC by using MPD for the analysis of two HCC tumor tissues that responded differently to galunisertib treatment.

Results
To investigate the molecular targets and druggable pathways for patients with HCC, we performed MACE genome-wide transcriptome sequencing and Exome-Seq all exome sequencing of tumoral (called 'tumor') and peritumoral (called 'normal') tissues extracted from tumors that had previously been shown to react differentially to TGFB1 and to the TGFBR1 inhibitor galunisertib. 15 Figure 1 Up-or downregulated genes in HCC tissues after galunisertib treatment. In (a), the total number of significantly (log 2FC 42 and P-value o0.01) up-or downregulated genes in responder-and non-responder tumor tissues, in comparison with corresponding normal tissues. In (b), the total number of overexpressed genes in both responder and non-responder tumor tissues compared with their corresponding normal tissues. The Venn diagram presents the total number of dissimilar and shared genes in responder and non-responder tumor tissues based on the total number of differential genes. Within parentheses, their percentages are given in relation to the aggregate number of differentially expressed genes in both tumor tissues. In (c), MACE profiles of responder-and non-responder normal and tumor tissues differentiating between all tissues. Unsupervised hierarchical cluster analysis of 164 strongly expressed genes (based on high intensities (50% 4100) and high variability (interquartile range (IQR)41.5)) clearly separates normal and tumor tissues. Normalized expression values are rescaled as shown in the sidebar, where a positive number (red) indicates high expression and a negative number (blue) low expression Using MACE we found a large number of genes differentially expressed between normal and tumor samples of both patients. In comparison with the non-responder tissue, there were relatively higher number of significantly (absolute log 2 fold change (log 2fc) 42 or o − 2 and P-value o0.01) differentially expressed genes in the responder sample, as depicted in Figure 1a. A total of 403 genes were similarly upregulated in both tumors in comparison with both normal tissues (Figure 1b). Unsupervised hierarchical clustering of strongly expressed genes shows a clear distinction between the transcriptomes of normal and tumor tissues (Figure 1c). Principal component 1 explained 79.8% of the total variance in the entire transcriptome data estimated by PCA analysis of 164 strongly expressed genes (data not shown).
Genes differentially expressed in both responder and nonresponder tumor samples in relation to their respective normal samples were subjected to Gene Ontology and Reactome pathway analyses to identify molecular processes and pathways differentiating normal from neoplastic tissues. The most significant pathways that clearly differentiate the two tumors from their respective reference tissues are shown in Figure 2.
In the responder tumor all four genes in the reactome pathway 'regulation of TP53 expression' including PRDM1 (known also as BLIMP1) were upregulated. Thus, upregulation of this pathway is involved in abnormal growth of these cancer cells. The transcriptome of the responder tumor also showed enriched pathways related to the formation, degradation and remodeling of the extracellular matrix (ECM) indicated by reactome pathways 'dissolution of fibrin clot', 'degradation of the extracellular matrix', 'activation of matrix metalloproteinases', 'collagen degradation' and 'non-integrin membrane-ECM interactions'. These pathways together included the largest number of upregulated genes and may be crucial for cell migration and evasion of cancer cells from the tumor and hence metastasis. In the non-responding tumor, enrichment of the pathways 'loss of function of SMAD4 in cancer' and 'SMAD4 MH2 domain mutants in cancer' indicates that the non-responder tumor shows impaired TGFB1-related transcription initiation, which corresponds to the observed (nonresponder) phenotype. More than 70% of the genes in these pathways are upregulated in the tumor. Further, the transcriptome of the non-responding tumor tissue showedlike that of the responder tumorsome enriched ECM-related pathways (e.g. 'collagen degradation', 'activation of matrix metalloproteinases'). Other pathways such as 'PERK regulates gene expression' and 'ATF4 activates genes' indicate marked membrane stress and unfolded protein responses. Cell-cycle control by P53 and EMI1 is indicated by the enrichment of the pathways 'TP53 regulates transcription of cell cycle genes' 'TP53 regulates transcription of genes involved in G2 cell cycle arrest' and 'phosphorylation of EMI1'. Thus, upregulation of the NIMA kinase pathway may be the upstream signal for the observed induction of the immune response-related pathways 'CLEC7A/inflammasome pathway' and 'interleukin-19 (IL-19), -20, -22 and -24' that indicate enhanced innate defense and tissue repair processes in the non-responder tumor.
As depicted in Figure 3a in the responder tissue, there are more than 50 upregulated kinases engaged in at least one upregulated reactome pathway. To investigate whether these kinases would interact with each other and to help with selecting, especially promising drug targets with large numbers of interactors, we performed a STRING analysis of protein-protein interactions of upregulated kinases in the responder tumor. This analysis delivered a kinase interaction network in which i. A. SRC, MAPK1, FGR, SYK and LYN form strongly interconnected hubs ( Figure 3b).
Since kinase activity is usually regulated by phosphorylation and not by overexpression, upregulation of a kinase hints at an even more important role in cancer development and maintenance. The MPD pipeline offers the options to select drugs according to general reactome profiles and/or to focus on particular pathways such as immune-related pathways that are already suggested by their molecular signatures. The first option interrogates the Human Protein Atlas (http://www. proteinatlas.org) database containing~670 FDA-approved drugs targeting functions enriched in upregulated genes and the respective proteins facilitating that function. The complete list resulting from relating the responders tumors' individual molecular patterns to this database is depicted in Table 1.
A list of kinases that are upregulated in the responder tissue are reported in Table 2. Several of them are targeted by approved or experimental tyrosine kinase inhibitors in different stages of development, as well as by natural, low-toxicity anticancer compounds.
Correlating drug reactivity profiles to enriched reactome pathways for the non-responder yielded the list of drugs depicted in Table 3. The first drug on the list is marimastat, which blocks several matrix metalloproteinases (MMPs) that are upregulated in the tumor tissue. Since MMPs are needed for the dissolution of the ECM enabling cancer cells to leave the tumor, their inhibition could help to prevent metastasis. The COX-2 inhibitor celecoxib, however, which was developed for treating rheumatic and degenerative diseases, also inhibits a set of at least six proteins that are upregulated in the tumor.
As depicted in Table 4, the number of kinases engaged in upregulated pathways in the non-responder tissue and resulting STRING protein interaction network (not shown) is much smaller than in the responder tumor. Consequently also the number of kinase inhibitors potentially available for treatment of the tumor is reduced. The small list of upregulated kinases in the non-responder tumor (Table 4) delivers only a few addressable targets. There are three approved drugs available for inhibition of MAP2K6 that is engaged in three enriched pathways and also MAPK11 that functions in two enriched pathways is targeted by Regorafenib. However, whereas the latter drug targets several kinases of expected relevance in the responder tumor, it only targets one (MAPK11) in the non-responder tumor tissue. Additional studies are needed to learn whether the drug would have similar effects in both patients.
The upregulation of different immune-related pathways in both tumors (Table 5) indicates a potential for tailored treatment alternatives or amendments to pathway inhibitors by addressing the immune checkpoints in the two tumors. Therefore, in Table 5, we present a list of immune-related pathways enriched in upregulated genes in the two tumors.
In the responder tumor the 'PD-1 signaling' pathway is most prominent, indicating that immune responses in the tumor are under the control of the immune checkpoint inhibitor PD1. By contrast, the non-responder tissue does not seem to rely on this type of control of the immune system since none of the PD1 signaling pathway genes are upregulated ( Table 5). The second most upregulated immune pathway in the responder tumor is 'regulation of innate immune response to cytosolic DNA', which corresponds to the concomitantly upregulated 'stimulator of interferon gene (STING)-mediated induction of host immune responses' pathway also dealing with the modulation of immune responses by cytosolic DNA. STING is a cytosolic receptor that senses both exogenous and endogenous cytosolic DNA. It activates TANK-binding kinase 1/interferon regulatory factor 3, nuclear factor-κB (NF-κB) and signal transducer and activator of transcription 6 (STAT6) signaling pathways to induce robust type I interferon and proinflammatory cytokine responses. The STING signaling pathway can be targeted by STINGVAX. STINGVAX are cyclic dinucleotides that are ligands formulated with granulocytemacrophage colony-stimulating factor-producing cellular cancer vaccines. 'cytotoxic T-lymphocyte-associated protein 4 (CTLA4) inhibitory signaling' is another pathway that is highly enriched in upregulated genes in the responder tissue. The CTLA-4 is the inhibitory checkpoint receptor on T cells. CTLA-4 upregulation by certain Tregs is linked to host immune tolerance after liver transplantation.
In the non-responder tumor the 'CLEC7A/inflammasome pathway' is the most prominent dysregulated pathway potentially available for therapeutic intervention (Table 5). CLEC7A (also called dectin-1) is located on the surface of particular immune cells. It belongs to a class of C-type lectin  Besides MACE transcription profiling we performed comparative all-exome sequencing of the corresponding tumor and normal tissues of the responder and non-responder tumors to disentangle mutations that may eventually help to explain the differential gene expression in the two tumors. Moreover, information about mutated genes might also serve to avoid administering drugs targeting non-functional or irreversibly activated kinases and phosphatases. The number of SNPs/InDels in the tumor tissues versus the respective peritumoral tissue according to their genomic location in relation to genes is summarized in Table 6.
To identify potential mutated hub genes within the tumors that might interact with other mutated genes we again performed a STRING protein-protein interaction network analysis of genes carrying protein-changing mutations, Figure 4a (responder) and B (non-responder).
Like the responder interactions, also the non-responder STRING protein interaction network shows the RANBP2 protein as a central hub in the network of mutated proteins. As in the responder, also in the non-responder it interacts with several mutated zinc-finger proteins (ZNF 507, 484, 75D, 667, 540, 292, 750A) (Figure 4), although these differ from those of the responder. In addition, the TOP2A gene, another hub gene in the network with which it interacts, is also mutated in the non-responder. RANBP2 sumoylates the TOP2A protein in mitosis, which is required for the proper localization of TOP2A to centromeres.
Our assumption that TGFB function is impaired in the nonresponder tumor is not only based on its phenotype but also on the fact that this tumor carries a nonsynonymous C to A mutation in the TGFBR3 gene (also known as betaglycan) on chromosome 1 at position 91 716 677 (P-value = 0.000980), which, to our knowledge, has not been previously described in HCCs. Contrary to TGFBR1 and TGFBR2, which bind directly  Indicated as an add-on or prophylactic oral medication in the chronic treatment of mild atopic asthmatic children. Also used as self-medication for the temporary relief of itching of the eye due to allergic conjunctivitis (ophthalmic).
Abbreviations: FDR, false discovery rate; i.v., intravenous; JRA, juvenile rheumatoid arthritis; MS, multiple sclerosis; NSAID, non-steroidal anti-inflammatory drug; OA, osteoarthritis; RA, rheumatoid arthritis; RR, relapsingremitting; SP, secondary progressive Personalizing cancer therapy in HCC R Agarwal et al to some TGFB superfamily ligands, TGFBR3 is a coreceptor not only for TGFB but also for related factors such as activins, inhibins, growth differentiation factors and bone morphogenetic proteins. To answer the question whether in our particular (non-responder) tumor TGFB3 may also interfere with EMT, we took a closer look into the expression of TGFBrelated transcripts in both tumors. In the responder tumor and non-responder tumor in which we list TGFB pathway transcripts that are significantly differentially expressed (log 2fc ⩾ 2), we found much more upregulated transcripts in Abbreviations: DTC, differentiated thyroid cancer; mCRC, metastatic colorectal cancer; N/A, not applicable; TKI, tyrosine kinase inhibitor TKIs targeting the respective kinases are also indicated  As shown for pancreatic cancer, suppression of EMT concomitant with galunisertib administration could increase the sensitivity of the the responder tumor to the drug and thus would be an interesting option for a combination therapy. 19

Discussion
Personalized therapy is one of the biggest challenges to overcome successfully the heterogeneity of HCC and to offer patients the most effective treatment for those recipients who could not receive LT, LRTs or hepatic resection. 20 This study firstly demonstrates that tailoring treatment strategies according to the individual HCCs genetic profile provides a wealth of therapeutic choices going far beyond standard sorafenib treatment for liver cancers. Here, we demonstrate that neither of the HCC tissues, responder and non-responder to galunisertib effectiveness in vitro according to our previous investigation, 15 contained the somatic mutations most frequently found in HCCs, such as TERT or TP53 mutations. 21 However, the mutation that we suspect has the largest impact on the differential character of molecular patterns in the two tumors is the mutation in the TGFBR3 gene in the nonresponder tumor. In HCCs as in other cancers such as ovarian cancer, TGFBR3 may function as a cancer suppressor, 12 and in pancreatic cancer loss of TGFBR3 expression promoted cancer progression. 22 In a pancreatic cancer model of the EMT, TGFBR3 suppresses the associated increased motility and invasiveness. Suppression of motility and invasiveness does not depend on its cytoplasmic domain or its coreceptor function but is mediated by ectodomain shedding, generating soluble TGFBR3. 23 Since it has been shown that AXL activates autocrine TGFB signaling in HCCs, resulting in a poor prognosis, it is tempting to speculate that in the nonresponder tumordue to the mutation in the TGFBR3 genethe TGFB autocrine loop including AXL is blocked. 24 This is consistent with our previous observation whereby TGF-β circulating levels did not correlate with the staging of the disease. 25 The observed mutation of the TGFBR3 gene in the non-responder tumor may be indicative of a certain type of HCC and at the same time may serve as a biomarker for companion diagnostics of TGFBR1-targeting drugs like galunisertib. 26 Moreover, these results indicate that TGFBR3 may be a similarly promising drug target (not necessarily inhibitor) as compared with the TGFB receptors, because its mutation seemingly has pleiotropic effects on TGFB-related pathways. The imbalance of the TGF pathway may evoke several downstream effects, including the proteolytic remodeling of the ECM proteins by MMPs and the EMT, both leading to the progression of the cancer. 27 The modulation of the immune response has been reported also to be under the control of TGFB, although here we cannot decide whether the marked differences between immune-related reactome pathways that also distinguish the responder from the nonresponder tumor are another consequence of the difference in TGFB signaling or not. In any case, the differential control of the immune system by the two tumors would have a major impact on a potential combination therapy if it should include immune-modulating drugs. These showed promising results and even a complete response in cancers with an otherwise poor prognosis. 27 Thus, the responder would probably profit most from blocking PD-1 checkpoint inhibition by drugs such as nivolumab, ipilimumab or pembrolizumab, if necessary supported by STINGVAX treatment. In lung cancers, anti-PD-1 therapy profited from increased nonsynonymous mutations, especially in DNA repair pathways in comparison with tumors with lower mutation frequencies. It was reported that 69% of lung cancer patients with a high mutation frequency experienced durable clinical benefits from PD-1 blockade as compared with 20% of patients with fewer mutations. 28 High mutation rates were also significantly associated with progression-free survival. Similar results were reported by Rizvi et al.,29 who noted that this favorable result was at least in one case related to neoantigen-specific CD8+ T-cell responses, suggesting that anti-PD-1 therapy enhances neoantigen-specific T-cell reactivity. The number of protein changing mutations for our HCC tumors was 275 in the responder and 226 in the non-responder tumor. Thus, both our tumors contain considerably more somatic mutations than the lung cancers studied by the above authors. This larger number of mutations opens up the possibility to load synthetic RNAs coding for a large number of different tumor-derived neoantigens via intravenously administered RNA-lipoplexes onto the immune system where they may evoke effector and memory T-cell responses, and mediate IFNα-dependent rejection of tumors. 30 Alternatively, synthetic naked DNA can be administered to skin-resident dendritic cells via micropinocytosis, which then prime antigen-specific CD8+ T cells. 31 Whereas the responder tumor engages PD-1 checkpoint inhibition, the non-responder tumor has already upregulated the 'CLEC7A/inflammasome pathway'. Thus, the nonresponder would probably be best treated using a combination therapy that boosts its already existing readiness to fight foreign invaders by such approaches as vaccination with the Bacille Calmette-Guerin tuberculosis, which in the case of this patient may not only help protect against pathogens but also against the tumor. 32 The upregulation of the CLEC7A pathway may be exploited for therapeutic, anticancer overactivation of the IL-6/STAT3 influenced innate immune system by components of pathogen cell walls such as fungal β-1,3 glucans or other immune-stimulatory substances. These would not only (over)activate the CLEC7A pathway but also specific, TLR cascades upregulated in the non-responder tumor. [33][34][35] This immune-stimulatory strategy could also benefit from the knowledge of the mutated proteins in the tumor. These could be used to synthesize candidate-mutated T-cell epitopes that may be identified using a major histocompatibility complexbinding algorithm for recognition by tumor-infiltrating lymphocytes. CLEC7A drives IL-1B biogenesis and maturation through a noncanonical caspase-8-dependent inflammasome in the host innate immune system. 36 The receptor dimerizes upon ligand binding and phosphorylation by kinases of the SRC family. Activation of CLEC7A finally leads to the activation of transcription factor NF-κB, which then induces the production of inflammatory cytokines and chemokines such as TNF, IL-23, IL-6 or IL-2. 37 Dectin-1 reduced hepatic fibrosis and hepatocarcinogenesis by negative regulation of TLR4 signaling pathways. 38 In conclusion, the study we present here has the aim of contributing to the discussion as to how therapeutic options arising from our growing ability to characterize individual cancers at the molecular level may best be used to the benefit of the individual patient.

Materials and Methods
Ex vivo HCC tissue profiling and DNA-RNA extraction. Tumoral and peritumoral tissues were isolated collecting the peritumoral at least at 4-5 cm by the limit of the lesion. Ex vivo assay was performed as described previously. 15 Briefly, human HCC samples were cultured for 48 h in serum-free condition in the presence of galunisertib and TGF-β. Tissues were stored in liquid nitrogen before RNA-DNA isolation. Total RNA and DNA were isolated using the AllPrep DNA/RNA Mini Kit (Qiagen, Hilden, Germany) according to the supplier's protocol.
Massive analysis of cDNA ends. MACE libraries were prepared from the two peritumoral and the two tumor tissues using the MACE Kit (GenXPro GmbH, Frankfurt Germany) according to the supplier's protocol as described. 16 In short, polyadenylated mRNA was extracted (Dynabeads mRNA Purification Kit; Life Technologies, Whaltham, MA, USA) from 5 μg total RNA (large RNA fraction) and reverse transcribed with biotinylated oligo (dT) primers. cDNA was prepared and fragmented to an average size of 250 bp by sonification using a Bioruptor (Diagenode, Seraing, Belgium). Biotinylated cDNA 3′ ends were captured by streptavidin beads and bound to TrueQuant DNA adapters (Waltham, MA, USA) provided in the kit. The libraries were amplified by PCR, purified by SPRI beads (Agencourt AMPure XP; Beckman Coulter, Brea, CA, USA) and sequenced (NextSeq 500; Illumina Inc., San Diego, CA, USA).
Bioinformatics analysis of MACE data. A total of~114.12 million MACE sequencing reads was obtained from the four cDNA sequencing libraries. All PCRduplicated reads identified by the TrueQuant technology were excluded from the raw data sets. The remaining reads were further quality trimmed and the poly (A)-tail was clipped off. Filtered reads were aligned to the human reference genome (hg38) using the bowtie2 mapping tool. 39 Unmapped reads were then aligned to the transcriptome assembly using the same mapping tool. Finally, the respective bam files of each MACE library were merged into respective single combined bam files using SAMtools. 40 The hg38 refSeq annotation GTF file, which includes coding genes as well as long noncoding RNAs (http://genome.ucsc.edu/cgi-bin/hgTables), was imported into the htseq-count annotation tool to annotate merged bam files. This yielded a set of combined enriched gene count data as output. Next, normalization of enriched gene count data followed by differential gene expression analysis between matched tumor samples and peritumor samples were performed to estimate expression levels of 41 636 genes/transcripts and 432 long noncoding RNAs, using the DeSeq R/Bioconductor package. To account for multiple testing, the false discovery rate (FDR) was estimated according to Klipper-Aurbach et al. 41 Genes with P-values o 0.01, and an absolute log 2 fold change |(log 2fc)| 4 2 were considered differentially expressed. Differentially expressed genes were uploaded to the Database for Annotation, Visualization and Integrated Discovery (DAVID) using an enrichment cutoff of FDR o 0.05. 42 Differentially expressed genes were further assigned to biological pathways using the Reactome database, a curated knowledgebase of biological pathways in humans. 43 Immune-related pathways were separated out from the entire list of reactome pathway terms based on evidence collected from recently published literature and terms were linked with immune regulation. Then, immune pathways in both responder and non-responder samples were ranked via estimating the percentage of affected genes in a particular pathway. The STRING database v.10 44 was used with standard parameters to discover interactions between proteins encoded by genes that were significantly up-or downregulated or somatically mutated in the tumor, respectively.
Identification of drugs targeting kinases in reactome pathways enriched in upregulated genes. Kinases engaged in pathways enriched with upregulated genes in tumor tissues were identified through mapping via their reactome ID, followed by extraction of kinase genes from mapped reactome IDs. Those kinases that were upregulated at log 2fc4+1 and with a P-valueo0.05 in both responder and non-responder tumor samples were chosen for further drug enrichment/ interaction. Drug interaction analysis was performed by interrogating the Drug-Gene Interaction Database (DGIdb), a database and web interface for ascertaining known and potential drug-gene relationships applying standard parameters. 45 Exome library preparation, sequencing and analysis. DNA was quantified using Qubit HS dsDNA assay (Life Technologies). Exomes were captured from a total of 50 ng of DNA using Illumina's Nextera Custom Target Enrichment Kit (Illumina) according to the Illumina Nextera Rapid Capture Enrichment Guide (August 2013) . A total of 2 ×~121.2 million 151 bp long paired-end (PE) reads of tumoral and peritumoral tissues was generated using a NextSeq platform (Illumina). Initially, raw reads were checked for adaptors and low-quality bases. The remaining 2 × 117.64 million PE reads were aligned to the human genome (hg38) reference sequence using the bwamem mapping tool with default settings. The Picard Toolkit (San Francisco, CA, USA) was used to preprocess mapped bam files in which SortSam, MarkDuplicate and BuildBamIndex commands were primarily used before calling for variants including SNPs and InDels. The standard GATK pipeline (http:// www.broadinstitute.org/gatk/) was used to call variants, in which the first step was the local realignment of reads around InDels and SNPs using IndelRealigner and RealignerTargetCreator commands, respectively, to curtail the number of mismatching bases across all the reads. Thereafter, a Mutect2 command 46 was used to identify somatic SNPs and the HaplotypeCaller command was used to detect somatic InDels. We removed low-quality SNPs by applying high-stringency criteria using the VariantFiltration command with filtering options (QDo2.0|| FS460.0||MQo40.0||HaplotypeScore413.0||MQRankSumo − 12.5||ReadPos-RankSumo − 8.0, where || indicates OR). In the same way, low-quality InDeLs were removed with filtering options (QDo2.0||FS4200.0||ReadPosRankSumo − 20). Remaining high-quality somatic SNPs and InDels were annotated using SnpEff 47 to assign them to different gene-related categories such as nonsynonymous, exonic and so on. The STRING database v.10 44 was used again with standard parameters to discover interactions between gene products (proteins) encoded by genes that harbor protein-changing SNPs/InDels.
Identification of potential therapeutic targets and corresponding drug candidates using the MPD drug prediction tool. Drug enrichment analysis was performed by calculating the significance of overlaps between a subset of 625 FDA-approved druggable genes (obtained from http://www. proteinatlas.org), which are significantly upregulated in one or the other tumor tissue (thresholds: P-valueo0.05, log 2fc41.0). Moreover, they are enriched in overall drug interactions for each compound as compared with all genes in the human genome. P-values denoting the significance of the enrichment were calculated using hypergeometric distribution. Drug interactions were annotated using data from DGIdb (http://dgidb.genome.wustl.edu).

Conflict of Interest
The authors declare no conflict of interest.