extfieldsofstudy: Medicine Chemistry

text:
Fisetin targets YB-1/RSK axis independent of its effect on ERK signaling: insights from in vitro and in vivo melanoma models

The anti-proliferative activity of dietary flavonoid fisetin has been validated in various cancer models. Establishing its precise mechanism of action has proved somewhat challenging given the multiplicity of its targets. We demonstrated that YB-1 promotes epithelial-to-mesenchymal transition and its inhibition suppressed tumor cell proliferation and invasion. The p90 ribosomal S6 kinase (RSK), an important ERK effector, activates YB-1 to drive melanoma growth. We found that fisetin treatment of monolayer/3-D melanoma cultures resulted in YB-1 dephosphorylation and reduced transcript levels. In parallel, fisetin suppressed mesenchymal markers and matrix-metalloproteinases in melanoma cells. Data from cell-free/cell-based systems indicated that fisetin inhibited RSK activity through binding to the kinase. Affinity studies for RSK isoforms evaluated stronger interaction for RSK2 than RSK1. Competition assays performed to monitor binding responses revealed that YB-1 and RSK2 do not compete, rather binding of fisetin to RSK2 promotes its binding to YB-1. Fisetin suppressed YB-1/RSK signaling independent of its effect on ERK, and reduced MDR1 levels. Comparable efficacy of fisetin and vemurafenib for inhibiting melanoma growth was noted albeit through divergent modulation of ERK. Our studies provide insight into additional modes of regulation through which fisetin interferes with melanoma growth underscoring its potential therapeutic efficacy in disease progression.

serine 380 residues are located in the linker region within the 'turn motif ' and the 'hydrophobic motif ' sequences of the kinase, respectively. The currently accepted model of RSK activation maintains that ERK1/2 activation results in the phosphorylation of threonine 573 in the CTKD of RSK. The activated CTKD then autophosphorylates RSK at the serine 380 residue. However, this site may also be phosphorylated by other kinases. In addition, ERK might also phosphorylate RSK at threonine 359 and serine 363 residues 5 . Alternatively, docking of PDKI at the phosphorylated hydrophobic motif phosphorylates serine 221 in the NTKD activation loop resulting in RSK activation 4,5 .
RSK2 was found to be an essential regulator in tumor promoter induced cell transformation 6 . Activated RSK2 protein levels are highly abundant in human skin cancer tissues compared with normal skin 7 . Studies show that RSK through differential regulation of pro-apoptotic protein Bad mediates a MAPK-dependent tumor-specific survival signal in melanoma cells 8 . Others have demonstrated that stimulated ERK pathway reduces the sensitivity of melanoma cell lines to cisplatin through activation of RSK1 9 . Expression profiling analysis revealed that ERK-activated RSK induces transcription of an effective pro-motile invasive gene program which results in modulation of extracellular and the intracellular motility apparatus. Thus RSK serves as a key effector, from which multiple highly coordinate transcription-dependent mechanisms originate for stimulation of cell motility and invasion 10 . Given the fact that inhibition of RSK isoforms has proven effective in blocking invasion and metastasis, RSK inhibitors are being investigated for their specificity and selectivity in different tumor types 11 . More recently, RSK inhibition was shown to block cell proliferation and protein synthesis in BRAF/MEK inhibitor-resistant melanomas, establishing this pathway as a viable therapeutic strategy against chemoresistance 12 .
The dietary flavonoid fisetin is gaining much attention for its pleiotropic effects in various biological systems 13 . Fisetin has been studied for its neuroprotective activity and was shown to confer protection in multiple models of Huntington's disease through ERK activation 14,15 . We and others have delineated its anti-cancer properties in several tumor models including melanoma and prostate 13 . We demonstrated that fisetin inhibits melanoma growth through suppression of Akt/mTOR pathways 16,17 . Analyzing the mechanism(s) of fisetin-induced cytotoxicity, we showed that fisetin mediated cell death is associated with induction of ER stress and activation of apoptotic pathways in human melanoma cells 18 . In other studies, fisetin was found to potentiate sorafenib-induced apoptosis and abrogate tumor growth in mice implanted with BRAF-mutant melanoma cells 19 . Fisetin suppressed human melanoma cell invasion by inhibiting epithelial-mesenchymal transition (EMT) and augmented anti-invasive and anti-metastatic potential of sorafenib in melanoma cells 20 . Recently, we identified YB-1, a DNA-RNA binding protein with a role in EMT and chemoresistance, as a key target of fisetin in prostate cancer cells 21 . Here we extended these studies and analyzed the mechanism through which fisetin exerts its effect on YB-1/RSK axis in human melanoma cells.
Cell culture/treatment. WM35 and A375 cells were obtained from ATCC (Manassas, VA), NCI/ADR-Res cells from NCI (Developmental Therapeutics Program, Frederick, MD) while 451Lu cells were provided by Dr. Meenhard Herlyn (Wistar Institute, PA). Authentication report was provided by ATCC and NCI for post-freeze viability, growth properties, morphology, mycoplasma contamination, species determination (cytochrome c oxidase I assay and STR analysis), sterility test and human pathogenic virus testing. Cells were frozen in aliquots in liquid nitrogen and cultured within six months. 451Lu cells were authenticated upon receipt using STR analysis and all cell lines are regularly tested for mycoplasma contamination using MycoAlert Mycoplasma Detection Kit from Lonza (Basel, Switzerland). WM35, A375 cells were cultured in DMEM, NCI/ADR-Res in RPMI and 451Lu in MEM with 10% FBS and 1% penicillin-streptomycin, at 37 °C with 5% CO 2 in a humid environment. After treatment of cells with/without fisetin and/or inhibitors, whole cell lysates were prepared and western blot analysis was performed as described elsewhere 22 . Viability studies. Cells seeded in 24-well plates were treated at 70% confluence with/without fisetin and/or vemurafenib for 24 hours. Cell viability was determined by MTT assay as described 17 . Immunoprecipitation. Briefly, specific antibody was added to 100 μg of cell lysates after pre-clearing with protein A/G-agarose beads. After overnight incubation, protein A/G agarose beads were added and left for 3 h. After centrifugation and several washes, the pellet was resuspended with 2X sample buffer and electrophoresed on a SDS-PAGE gel.
Immunochemistry. For immunohistochemistry, sections from melanoma constructs, fixed in 10% formalin were deparaffinized and rehydrated in xylene and a graded series of ethanol and blocked in 2% goat serum/ PBS for 30 min. After incubation with specific primary antibody overnight and biotinylated secondary antibody for 1 h, slides were counterstained with hematoxylin and immunoreactive complexes were detected using 3,3′-diaminobenzidene (Dako Corp., CA). Immunocytochemical analysis was performed after seeding cells in 2-chamber tissue culture slides treated with/without fisetin and probing with fluorophore linked secondary antibodies 18  Kinase binding assay. The KINOMEscan screening platform (San Diego, CA) was employed to quantitatively measure interactions between fisetin and RSK isoforms. Binding constants (Kds) were calculated by measuring the amount of kinase captured on solid support as a function of the test compound concentration, using the Hill equation. The Hill Slope was set to −1. Curves were fitted using a non-linear least square fit with the Levenberg-Marquardt algorithm 17 . Kinase activity assay. In vitro kinase activity profiling of RSK isoforms was performed at Reaction Biology Corporation (Malvern, PA). Briefly, the specific kinase/substrate pair alongwith required cofactors was prepared in reaction buffer. Fisetin/staurosporine dissolved in DMSO was added to the mixture, followed by 32 P-ATP. After incubation for 2 h, reactions were spotted onto filter paper which binds the radioisotope-labeled catalytic product 17 . After subtraction of background derived from control reactions containing inactive enzyme, kinase activity data was expressed as % remaining kinase activity. Binding of analyte to the immobilized protein was monitored in real time to obtain on (ka) and off (kd) rates. The equilibrium constant (K D ) was calculated from the observed ka and kd or by steady state kinetics as for the small molecules. Stocks were prepared in 100% DMSO, and further dilutions were made in assay buffer containing 10 mM HEPES buffer (pH 7.4), 150 mM NaCl, 3 mM EDTA, 0.05% P20 (polyoxyethylenesorbitan) and 0.1% BSA. Full kinetics analysis was performed at using analyte concentrations 50 nM (or 100 nM without BSA) to 0, run with serial dilutions, 50, 25, 12.5, 6.25, 3.125, 0 nM, and a flow rate of 30 ul/min. For competition assays, fisetin alone (1 µM and 10 µM) and RSK2 alone (100 nM) were flown over immobilized YB-1 to monitor binding response. The concentration of RSK2 was kept constant, but pre-incubated with increasing amounts of fisetin and then injected to determine the effect of fisetin on binding activity of RSK2 to YB-1. An aliquot of cell lysate was spiked with PLX4720 (25 ng) [internal standard] and alkalinized by addition of two volumes of pH 11 buffer (1 mM sodium hydroxide and 0.5 mM sodium bicarbonate). Organic layer was separated and dried under nitrogen. Samples were reconstituted in 100 μl of mobile phase and transferred into HPLC glass vials. Chromatographic analysis was performed using the AQUITY UPLC system (Waters, Milford, MA). Chromatographic separation was achieved using an Agilent Technologies Eclipse XDB-C18 column (4.6 × 50 mm) with 1.8-μm Zobrax Rx-SIL as the stationary phase. Mobile phase consisted of 20 mM ammonium formate with 0.1% formic acid and acetonitrile (30:70 v/v), delivered at a flow rate of 0.25 ml/min. Samples were analyzed using an electrospray probe in the negative ionization mode operating at a voltage of 2.96 kV for both vemurafenib and PLX4720. The assay was sensitive and linear over a range of 1.2 ng/ml-1.2 μg/ml, with the coefficient of variation being less than 15% over the entire range.
Xenograft studies. Athymic (nu/nu) female nude mice (Harlan, Madison, WI), 4-5 weeks old, were housed under pathogen-free conditions, in the University of Wisconsin Animal Resource Facility, with a 12-hour light/12-hour dark schedule and fed an autoclaved diet ad libitum. All animal experiments were approved by and carried out according to Institutional Animal Care and Use Committee guidelines, University of Wisconsin-Madison, protocol # M005212. Briefly, 1 × 10 6 451Lu cells were injected s/c into each flank of 18 mice randomly divided into 3 groups. The first group received intraperitoneal injection of fisetin 40 mg/kg b.wt.
(1 mg/animal) thrice weekly, the second group received vemurafenib by oral gavage 80 mg/kg b.wt. (2 mg/animal) dissolved in DMSO, 5 days/week, while the third group served as the DMSO treated control group. Mice weight and tumor volumes were recorded twice weekly. The animals were sacrificed when tumors reached a volume of 1200 mm 3 in the control group. Samples were collected and stored at −80 °C until further analysis.

Fisetin inhibits YB-1 in BRAF mutant melanoma cells. YB-1 expression is reportedly upregulated
during melanoma progression and has been linked to proliferation, invasion and resistance to apoptosis 23 . We examined the effect of fisetin on endogenous levels of YB-1 in highly aggressive BRAF mutant melanoma cells. Fisetin treatment (0-80 µM:24 h) to A375 and 451Lu melanoma cells resulted in a dose-dependent decrease in YB-1 phosphorylation concomitant with downregulation of total protein levels, as assessed by western blot and immunocytochemical analysis (Fig. 1A,B, Sup. Fig. 1A-C). We selected WM35 melanoma cells, representative of non-invasive radial growth phenotype with an intrinsically low YB-1 expression for overexpression studies (Sup. Fig. 1D). A similar effect was observed in WM35 melanoma cells where exogenous YB-1 levels were significantly decreased with fisetin treatment (Fig. 1E). We then examined the effect of fisetin on YB-1 mRNA levels. qPCR analysis showed a significant reduction in YB-1 transcript levels in fisetin-treated A375 cells (Fig. 1D). Finally, we studied the expression of YB-1 in 3-D melanoma constructs at different stages of melanoma progression. Immunohistochemical analysis of untreated melanoma constructs showed increased expression of YB-1 which correlated with melanoma progression. In contrast, fisetin-treated melanoma constructs, at both time points (days 12 and 16) showed diminished staining of YB-1 versus the untreated controls (Fig. 1E). We had previously shown that forced expression of YB-1 promoted EMT in prostate cancer that was effectively inhibited by fisetin 21 . We examined EMT markers in fisetin-treated A375 melanoma cells and observed a decrease in vimentin and N-cadherin accompanied by an increase in E-cadherin. Moreover, fisetin treatment resulted in downregulation of MMP-2 and MMP-9 (data not shown) in monolayer and 3-D melanoma cultures (Sup. Fig. 1F,G).
Fisetin binds to RSK and suppresses its kinase activity. RSK has been recognized as a novel activator of YB-1, able to phosphorylate its serine 102 residue 24 . To examine if suppression of YB-1 phosphorylation by fisetin is linked to its inhibitory effect on RSK, we utilized radiometric assays (Reaction Biology Corp) to measure the kinase catalytic activity of RSK isoforms. Comparative studies in the presence of 32 P-γ-ATP showed that fisetin inhibited RSK3>RSK2>RSK1 with IC 50 values ranging from 8.91E-07, 2.78E-06 and 3.79E-06 M, respectively ( Fig. 2A, Sup. Fig. 2A). Employing the Kinomescan screening platform, ATP site-dependent competition binding assays were performed to quantitatively measure interactions between fisetin and RSK. Screening hits were identified by measuring the amount of kinase captured in treated versus control samples using qPCR, subsequent to which dissociation constants (Kds) for compound-kinase interactions were calculated. Remarkably, fisetin exhibited the highest affinity for RSK2 at its NTKD (50 nM) while the Kds for RSK1 and RSK3 ranged from 8100:2500 and 5200:5800 nM for N and C-termini, respectively (Fig. 2B,C, Sup. Fig. 2B). Fisetin/RSK2 complex augments its binding to YB-1. Our in silico studies indicated that fisetin binds to YB-1 on residues β1-β4 strands within its cold shock domain 21 . To examine if binding of fisetin to RSK isoforms had any impact on its binding to YB-1, we performed real-time interaction analysis using SPR assays. Binding evaluation was conducted by immobilizing the YB-1 protein on a sensor surface. The analytes (RSK1/RSK2) were then injected in solution over the surface. Changes in SPR response were analyzed, providing kinetic and affinity data. The scouting full kinetics performed to evaluate binding between YB-1 and RSK1/YB-1 and RSK2 indicated stronger interactions for the latter (K D = 2.46E-04M for RSK1 vs 5.34E-07M for RSK2, respectively) (Fig. 3A). Subsequent full kinetic assays repeated for RSK2/YB-1 in the presence of BSA, yielded K D values very close to those without BSA (4.02 × −07M) validating the system (Fig. 3B). Competition assays were carried out flowing first fisetin alone and then RSK2 alone over immobilized YB-1 to monitor the binding response (data not shown). In the next step, RSK2 and fisetin were preincubated, and binding experiments were finalized (Fig. 3C). Contrary to our expectations, data obtained revealed that RSK2 and YB-1 do not compete with each other; rather binding of fisetin to RSK2 seems to promote its binding to YB-1. Reverting to a cell-based system, we treated A375 melanoma cells with RSK inhibitor BI-D1870 with or without fisetin and scrutinized YB-1 expression. Both fisetin and BI-D1870 had similar effect on the phosphorylation of RSK at 359/363 residues. However, in comparison to fisetin BI-D1870 had only a modest effect on YB-1 expression (Fig. 3D, Sup. Fig. 3). This is in agreement with previous observations where BI-D1870 decreased phosphorylation of YB-1 but total protein levels remained unaltered 25 . In our studies, the combinatorial treatment did not yield an additive decrease in YB-1 levels. Taken together, the data suggest that formation of YB-1/RSK complex may be an important regulatory event in fisetin-mediated decrease of YB-1.  . Briefly, YB-1 was directly immobilized on the sensor chip by amine coupling and RSK1 or RSK2 was flown over the proteincoated chip at different concentrations (100-0 nM). ka and kd or steady state kinetics were used for determining the K D . Chi square (χ 2 ) analysis was carried out between the actual sensorgram (colored line, bottom) and the sensorgram generated from the BIAnalysis software (black line, top) to determine the accuracy of the analysis. χ 2 value within 1-2 was considered significant (accurate), and below 1 was highly significant (highly accurate).   26 . Moreover, increase in MDR1 expression has been linked to YB-1 phosphorylation and nuclear translocation in various tumor models 27 . Our studies demonstrated significantly diminished MDR1 levels in A375 and 451Lu melanoma cells (Fig. 4A, Sup. Fig. 4A). The NCI/ADR-Res is an ovarian tumor cell line that expresses high levels of MDR1. We therefore examined the susceptibility of NCI/ADR-Res cells to fisetin. Consistent with our previous data, these cells showed reduced viability and decrease in MDR1 and YB-1 levels upon treatment with fisetin ( Fig. 4B,C, Sup. Fig. 4B-D). Time course analysis in A375 cells revealed decreased MDR1 protein expression as early as 12 h post treatment. A parallel decrease in RSK phosphorylation at threonine 359 /serine 363 residues was observed (Fig. 4D, Sup. Fig. 4E). Evaluation of fisetin-treated 451Lu and NCI/ADR-Res cell lysates, incubated with agarose beads coated with antibody against RSK validated the presence of YB-1/RSK complex and interaction between RSK and MDR1 (Fig. 4E, Sup. Fig. 4F). Immunocytochemical analysis demonstrated enhanced staining of YB-1 and MDR1 in untreated melanoma cells that was significantly reduced with fisetin treatment (Fig. 4F). To further study these interactions, we overexpressed YB-1 in WM35 melanoma cells and assessed MDR1 expression. Interestingly, YB-1 overexpression increased MDR1 expression yet had little effect on phosphorylated RSK levels suggestive of a hierarchy where YB-1 functions downstream of RSK. Regardless, fisetin treatment abrogated RSK phosphorylation at threonine 359 /serine 363 and decreased MDR1 levels in YB-1 overexpressing WM35 melanoma cells (Fig. 4G, Sup. Fig. 4G).
We hypothesized that the strong suppressive effect of fisetin on MDR1 expression should potentially translate into potent inhibition of MDR1 substrate efflux, enhanced intracellular accumulation and therapeutic efficacy. To corroborate this line of thought, we tested whether fisetin can augment the inhibitory effect of BRAF inhibitor vemurafenib (PLX4032), a known MDR1 substrate. After establishing the IC 50 of vemurafenib (Sup. Fig. 4H) and fisetin (data not shown) in melanoma cell lines, we evaluated the combinatorial effect of vemurafenib and fisetin on melanoma cell viability. Our studies showed that combination of fisetin and vemurafenib was significantly more potent in decreasing cell viability in both A375 (data not shown) and 451Lu melanoma cells when compared to treatment with each compound alone (Fig. 4H). We then asked if the increased efficacy of the combination was due to increased cellular accumulation of vemurafenib, due to the inhibitory effect of fisetin on MDR1 mediated drug efflux. To answer this question, 451Lu (data not shown) and A375 melanoma cells, harvested 4 h and 8 h post treatment with fisetin and vemurafenib, were analyzed by LC/MS/MS (Fig. 4I). The combination group compared to cells treated with vemurafenib alone did not have any appreciable increase in vemurafenib concentration. To further study the functional relevance of fisetin-mediated MDR1 inhibition, we evaluated the efflux of MDR1 substrates DiOC2(3) and rhodamine in A375 melanoma and NCI/ADR-Res ovarian cancer cells under physiologic conditions. The relative fluorescence of cell populations treated with fisetin or known MDR1 inhibitors was quantified. Increased fluorescence in cells treated with MDR1 inhibitor vinblastine, indicated effective blockage of DiOC2(3) and rhodamine (data not shown) efflux. Fisetin on the other hand had only a modest effect on blocking the efflux, and the increase in fluorescence was not significant when compared to DMSO treated cells (Fig. 4J). This indicated that fisetin is not effective in inhibiting the efflux of MDR1 substrates, despite its suppressive effect on RSK/YB-1/MDR1 signaling axis.

Fisetin mediated decrease in RSK/YB-1 does not require suppression of ERK signaling.
Given the regulatory role of ERK signaling on RSK/YB-1 axis and its importance in BRAF mutated melanoma, we examined the effect of fisetin on ERK phosphorylation. Intriguingly, 24 h post treatment, fisetin increased the phosphorylation of ERK1/2 at threonine 202 /tyrosine 204 residues in 451Lu (data not shown) and A375 melanoma cells, at doses linked to decreased cell viability (Fig. 5A, Sup. Fig. 5A). Moreover, fisetin treated YB-1-overexpressing WM35 melanoma cells (Fig. 5B, Sup. Fig. 5B) and chemoresistant NCI/ADR-Res ovarian cancer cells (Fig. 5C, Sup. Fig. 5C) also exhibited a dose dependent increase in ERK phosphorylation at 24 h. We further observed that fisetin mediated dose-dependent decrease in RSK phosphorylation at threonine 359 /serine 363 residues (Fig. 5A,C,D) was accompanied with increased phosphorylation at the serine 380 residue (Fig. 5A,C, Sup. Fig. 5A-C). To investigate if RSK phosphorylation was regulated by ERK, A375 melanoma cells were exposed to FR180204, a selective ERK1/2 inhibitor. A significant decrease in RSK phosphorylation at threonine 359 /serine 363 residues was observed when cells were treated with fisetin or FR180204. Inhibition of ERK did not interfere with RSK phosphorylation at threonine 359 /serine 363 residues in fisetin-treated cells (Fig. 5E,F, Sup. Fig. 5D,E). However, fisetin-mediated increase in serine 380 phosphorylation was abrogated with ERK inhibition suggesting that ERK was driving phosphorylation at this residue (Fig. 5E, Sup. Fig. 5D,E). Importantly, fisetin mediated decrease in YB-1 levels was unaffected by ERK inhibition (Fig. 5E,F). Thus, dephosporylation of RSK at threonine 359 /serine 363 residues seems to be an important determinant in fisetin mediated inhibition of YB1, independent of ERK signaling.

Comparable inhibition of melanoma tumor growth by fisetin and vemurafenib albeit through divergent regulation of ERK signaling. The biochemical affinity of vemurafenib for mutated BRAF
translates into potent suppression of ERK phosphorylation and inhibition of cell proliferation 28 . We compared the effect of fisetin and vemurafenib on RSK/YB1 signaling in in vitro melanoma cultures. A375 cells treated with vemurafenib for 24 h showed a comparable decrease in YB-1 expression and RSK phosphorylation levels as witnessed with fisetin treatment. In contrast to vemurafenib, fisetin treatment induced ERK1/2 phosphorylation in A375 melanoma cells (Fig. 6A, Sup. Fig. 6A). However, cell free kinase activity studies indicated that fisetin mediated inhibition of ERK1/2 kinase activity was greater than staurosporine, used as a positive control (Fig. 6B,  Sup. Fig. 6B). To evaluate this apparent discrepancy we performed a time course analysis of fisetin-treated A375 cells. We observed that an initial increase in ERK1/2 phosphorylation was followed by a sustained decrease, suggestive of a biphasic response (Fig. 6C, Sup. Fig. 6C). Next studies were conducted to ascertain if increase in ERK phosphorylation was driving fisetin-induced apoptosis in melanoma cells. We found that ERK inhibition did not abrogate caspase-3 cleavage in fisetin treated cells (Fig. 6D, Sup. Fig. 6D). We further compared the effects of fisetin and vemurafenib in athymic nude mice implanted with 451Lu melanoma cells. The animals were divided into three cohorts, each with 6 animals that were administered fisetin/DMSO (intraperitoneally) or vemurafenib (orally). A smaller average tumor volume was observed in mice treated with fisetin and vemurafenib. Upon termination of the study, tumor volumes evaluated between the three groups showed no statistical difference in tumor inhibition between animals receiving fisetin or vemurafenib (p = 0.2). Notably, a significant decrease in tumor volumes was observed in animals receiving fisetin or vemurafenib when compared to DMSO treated controls (Fig. 6E). In the control group, the average tumor volume of 1179.7 mm 3 was reached at day 35, while mice receiving fisetin (1 mg/animal) and vemurafenib (2 mg/animal) had an average tumor volume of 789.3 (p = 0.05) and 641.7 mm 3 (p = 0.002) respectively. Similar to our in vitro studies in A375 melanoma cells, 451Lu xenograft lysates treated with fisetin and vemurafenib showed dephosphorylated RSK and decreased YB-1 levels (Fig. 6F,  Sup. Fig. 6E). However, in contrast to our in vitro data, we did not observe increased ERK phosphorylation in fisetin-treated mice. Animals treated with fisetin showed reduced phosphorylation of ERK1/2 as was observed in 3-D melanoma cultures exposed to fisetin for 16 days (Fig. 6F, Sup. Fig. 6F).

Discussion
New treatment options are constantly being explored to improve the prognosis of malignant melanoma, an aggressive and notoriously chemoresistant cancer 29 . The specificity paradigm that ascribed central relevance to achieving target specificity of drug candidates is now being reconsidered. Indeed there is considerable evidence that polypharmacological drug behavior is often responsible for improved therapeutic efficacy. Studies show that active compounds may interact with multiple targets, which in most cases, are related to each other 30 . Our lab and others have investigated and established the efficacy of dietary flavonoid fisetin in inhibiting melanoma growth in various experimental models, both in vitro and in vivo 13 . Defining the precise mechanism of action of fisetin remains a challenge due to multiple targets that may reside within sequentially regulated phosphorylation cascades with several interacting partners. Ongoing studies have delineated various signaling pathways involved in growth and proliferation that are evident targets of fisetin 13 . We demonstrated previously that fisetin inhibited YB-1 and markers of EMT in prostate cancer cells. Here we validated the effect of fisetin on YB-1 and evaluated its interaction with RSK in BRAF-mutant melanoma cells. The insights gained from this study include (i) fisetin mediated inhibition of YB-1 extends to other cancer cell types including melanoma and ovarian cancer cells (ii) fisetin binds to RSK2 avidly and this complex augments the binding of fisetin to YB-1 (iii) fisetin suppresses MDR1 downstream of YB-1 but inhibition does not translate into decreased substrate efflux (iv) fisetin mediated RSK suppression is independent of its effect on ERK signaling (Fig. 7).
Disruption of YB-1 phosphorylation by site-directed mutagenesis suppressed tumor growth in various tumor models 31 . Cell permeable interference peptides acting as molecular decoys and kinase inhibitors are being explored to block activation of endogenous YB-1 32 . While a considerable proportion of its oncogenic effect is thought to be facilitated by its transcriptional activity mediated by serine 102 phosphorylation, YB-1 can suppress or stimulate translation. It was noted that enhanced cytoplasmic expression of YB-1 in breast cancer cells may result in induction of tumor dormancy which can be reverted upon re-phosphorylation of YB-1 33 . Thus targeting total YB-1 levels through gene silencing studies was suggested to be a better approach for decreasing YB-1 levels and inhibiting tumor growth and invasion 34 . Our studies in monolayer and 3-D melanoma cell cultures showed that fisetin treatment to melanoma cells achieves the desired decrease in YB-1 abundance by targeting YB-1 at both protein and mRNA levels (Fig. 1).
There exists significant cross talk between YB-1 and RSK, crucial signaling molecules implicated in tumor growth, EMT and invasion [35][36][37] . Small-molecules that modulate interactions between specific proteins represent a promising avenue for therapeutic intervention in a variety of settings. Binding of the compound to the proteins may result in change in conformation that is distinct from the unbound or protein-bound conformations 38 . In accordance with cell-free data that provided evidence that fisetin binds avidly to RSK and inhibits its kinase activity, we observed a decrease in RSK phosphorylation at threonine 359 /serine 363 residues in cells treated with fisetin. In fact, both fisetin and the RSK inhibitor BI-D1870 exerted analogous effects on RSK phosphorylation 25 . Future studies will describe whether fisetin and BI-D1870, an ATP-competitive inhibitor, share common binding sites or target different areas on the kinase. Employing SPR studies, we evaluated the binding of fisetin with different isoforms of RSK and found strong binding with RSK2. We observed that addition of fisetin in the presence of RSK2 facilitated its binding to YB-1. We postulate formation of a 'ternary complex' between RSK2/YB-1/fisetin, where fisetin stabilizes the interaction between RSK2 and YB-1 by acting at the interface between the two kinases. The interaction between RSK and YB-1 reportedly culminates in phosphorylation and activation of YB-1 24,39 . In our model, however, binding of RSK2 with YB-1 does not translate into increased activity of YB-1. NMR studies elucidating folding dynamics of proteins indicate that binding of a ligand at a specific site may bias the conformational equilibrium toward a structure in which the active site may or may not be active 40 . Thus, binding of the allosteric effector by altering the conformation of the protein may well result into an altered functional output 40 . It is tempting to speculate that fisetin in a complex with YB-1 and RSK induced allosteric modifications in YB-1 such that RSK is unable to phosphorylate it at the serine 102 residue. This hypothesis is consistent with our previous observations where we showed that fisetin not only binds within the CSD of YB-1 but also alters its topology 21 . At this stage we have not explored concomitant interactions, yet we cannot exclude that this complex may also interfere with phosphorylation of RSK at threonine 359 /serine 363 by upstream kinases. The fact that ERK mediated phosphorylation of RSK at serine 380 residue remained unaffected by this complex indicates that this effect may be specific to the afore-mentioned site. It is interesting to note that fisetin did not affect RSK protein levels as opposed to YB-1 where a reduction in phosphorylation was accompanied by suppressed protein and mRNA levels. This suggests that fisetin may regulate YB-1 not only through modulation of phosphorylation through Figure 7. Schematic diagram depicting the proposed mechanism of action of fisetin: SPR studies showed that interaction of fisetin with RSK2 promotes its binding to YB-1. Fisetin mediated downregulation of YB-1/RSK signaling is associated with a decrease in MDR1. Fisetin-treated cells showed an initial increase in ERK1/2 phosphorylation at threonine 202 /tyrosine 204 followed by a sustained decrease, suggestive of a biphasic response. Fisetin mediated decrease in RSK/YB-1 is independent of ERK. Inhibition of ERK did not interfere with RSK phosphorylation at threonine 359 /serine 363 residues in fisetin-treated cells. However, fisetin mediated increase in RSK serine 380 phosphorylation was abrogated with ERK inhibition suggesting that ERK was driving phosphorylation at this residue.
The well-characterized multidrug transporter MDR1 is a major component of the energy-dependent efflux system, responsible for removal of chemotherapeutics out of tumor cells 41 . Although indisputable that MDR1 effluxes xenobiotics from cells, insights into the biological function of MDR1 and related transporter molecules suggest additional roles in chloride channel activity, cholesterol metabolism and immune cell function. Its involvement in regulation of cell differentiation, proliferation, and survival has been documented 42 . MDR1 was found to protect cells from apoptosis induced by a wide array of cell death stimuli that relied on activation of intracellular caspases for full function 42,43 . Recently, it was shown that MDR1 associated with the calcium-dependent phospholipid-binding protein Annexin A2 to promote invasion in multidrug-resistant breast cancer cells 44,45 . In our studies we observed that fisetin downregulated MDR1 both at transcription and translational levels however the decrease in MDR1 was not associated with intracellular accumulation of MDR1 substrates. It is probable that MDR1 partakes in YB-1/RSK driven EMT and cell survival signals. Fisetin mediated suppression of melanoma cell proliferation and invasion may be linked to its inhibitory effect on YB-1/RSK signaling that in turn regulates MDR1 thereby inhibiting tumor growth.
An important aspect of our work was to elucidate the regulation of ERK signaling by fisetin in melanoma cells. Studies in 3-D and animal models representative of long term treatment response showed dephosphorylation of ERK with fisetin treatment. This correlated with data acquired from cell-free studies where fisetin was found to be a potent inhibitor of ERK activity. Interestingly, time course studies in monolayer melanoma cells indicate a biphasic response where suppressed ERK activity was preceded by an initial increase in ERK phosphorylation. We have previously demonstrated that fisetin treatment results in increased nitric oxide levels in melanoma cells 18 . In certain cell types, nitric oxide has been shown to promote ERK signaling via redox reactions 46 . It can be surmised that the initial increase in ERK observed in fisetin treated cells may be associated with enhanced nitric oxide levels, however additional studies are warranted. Importantly, our in vitro data unambiguously shows that suppression of RSK phosphorylation by fisetin is not mediated through negative regulation of ERK signaling. Moreover, fisetin-induced cell death is independent of ERK signaling as shown in Fig. 6. Comparative studies of fisetin with vemurafenib underline the efficacy of fisetin and shed light on the molecular mechanisms involved in inhibition of melanoma growth. Going forward a combinatorial study employing low doses of vemurafenib and fisetin can provide powerful evidence to examine the dietary flavonoid in pre-clinical trials and subsequently in clinical settings.


Exercise, Telomeres, and Cancer: “The Exercise-Telomere Hypothesis”

Telomeres are genomic complex at the end of chromosomes that protects the DNA and telomere length (TL) is related to several age-related diseases, lifespan, and cancer. On the other hand, cancer is a multifactorial disease that is responsible for reduce the quality of life and kills millions of people every year. Both, shorter TL and cancer are related and could be treated or prevented depending of the lifestyle. In this review we discuss the possible role of exercise in the relationship between shorter telomeres, telomerase activity, and cancer. In summary, there is evidence that exercise leads to less telomere attrition and exercise also may diminish the risk of cancer, these two outcomes are possible intermediated by a reduction in oxidative stress, and chronic inflammation. Although, there is evidence that shorter TL are associated with cancer, the possible mechanisms that one may lead to the other remains to be clarified. We assume that humans under cancer treatment may suffer a great decrease in quality of life, which may increase sedentary behavior and lead to increased telomere attrition. And those humans with already shorter TL likely lived under a poor lifestyle and might have an increased risk to have cancer.

INTRODUCTION
Telomeres are specialized structures at the ends of linear chromosomes (Kipling, 1995;Blackburn et al., 2015). A telomere is a region of repetitive nucleotide (organic molecule) sequences at each termini of a chromosome (DNA molecule), which protects the end of the chromosome to maintain genomic stability, avoiding degradation and fusion (Blackburn, 2010). Telomere's name is derived from the Greek nouns telos (τέλoς), that means "end" and meros (µέρoς, root: µερ-), that means "part." As a physiological cellular process, a small part of telomeric DNA is lost with each cell division (Shammas, 2011). Therefore, telomeres length (TL) shortens every time a cell division occurs, which makes it a marker of biological aging (Blackburn et al., 2015). TL is also associated with large number age-related disorders, such as type 2 diabetes, hypertension, Alzheimer's and Parkinson's disease, and even cancer (Shammas, 2011).
When TL is short enough, damage to the DNA may cause the cell to produce non-functional proteins, which may reverberate in two different paths: (i) the non-functional cell is led to apoptosis, which is roughly the senescence process that humans age; or (ii) the non-functional cell continues to produce non-functional proteins, which may become non-functional cell and eventually becomes a cancer (d'Adda di Fagagna, 2008;Maciejowski and de Lange, 2017;Cleal et al., 2018).
On the other hand, life-long exercise practice has being shown to reduce biochemical factors related to telomere attrition, such as oxidative stress  and chronic inflammation (Gleeson et al., 2011). There are also some studies reporting that a sedentary behavior is associated with shorter TL (Sjogren et al., 2014), and furthermore, elite athletes, and master athletes have longer TL than their non-athlete counterparts (Muniesa et al., 2017;Simoes et al., 2017).
Although literature is almost consensual pointing the benefits of exercise to prevent and treat a great variety of cancers, the relation between telomere-exercise-cancer still needs clarification. In this narrative review, we aimed to identify the beneficial effects of exercise to prevent and treat cancer. Moreover, we also aimed to identify which methodological aspects are missing in the literature that could help in the development of future experimental designs.

TELOMERE AND TELOMERASE DYNAMICS
Telomeres are specialized functional complexes (parts of DNA molecules) that are cited at the end of DNA; they protect the ends of eukaryotic chromosomes from terminal attrition (Blackburn, 2001). On the other hand, telomerase is an enzyme responsible for synthesizes DNA at the chromosome terminal, lengthening the terminal regions of telomeric DNA by RNA-templated addition of tandemly repeated (Blackburn, 2001).
Telomerase is an RNA-dependent DNA polymerase (ribonucleoprotein) that synthesizes telomeric DNA sequences and provides the molecular basis for unlimited proliferative potential (Cong et al., 2002). Furthermore it counteracts telomere erosion and under certain circumstances may also elongate TL (Aviv, 2002). Telomerases include a protein subgroup of specialized reverse transcriptase enzymes known as Telomerase Reverse Transcriptase (TERT). The modulation of telomerase activity may have significant anti-aging implications as well as anticancer therapy (Cong et al., 2002). A study previously reported that nucleotide polymorphisms in the TERT, are strongly associated with decreased breast cancer risk (Helbig et al., 2017).

TELOMERES AND AGING
Telomeres are the parts of chromosomes that affect aging, and their length (TL) appears to be an indicator of biological age (Sanders and Newman, 2013;Zhang et al., 2014;Jylhava et al., 2017). Telomeres become shorter due to cell division and lifestyle factors; critically shortened telomeres are linked to cellular dysfunction and aging (Garcia-Martin et al., 2017). Telomeres shorten with age, with each round of mitosis (division of the nucleus), due to the inability of the DNA replication mechanisms to read and copy the ends of linear chromosomes (phenomenon known as the "end replication problem"), and also from oxidative stress, which induces DNA damage (Eisenberg, 2011).
Telomeres act as a buffer against the loss of coding or regulatory DNA during cell proliferation (Helbig et al., 2017). Regarding the role of telomerase activity on cell aging, it has being shown that cells with a perfect regulated telomerase activity are considered immortal because they can divide past the Hayflick limit without entering in phase of senescence (biological aging) or apoptosis.
Shorter TL is also closely associated to almost all agerelated diseases (Blackburn et al., 2015), and the literature has been slashing to point that the main reason modulating a possible acceleration or attenuation in this process is the lifestyle (Blackburn et al., 2015;Denham et al., 2016). An inadequate lifestyle and aging are both closely related to chronic inflammation (Gleeson et al., 2011), increased oxidative stress  and decrease in telomerase activity (Blackburn et al., 2015).

TELOMERES AND EXERCISE
Lifestyle is a multifactor aspect with regards to general health and could be simplified as nutrition habits, stress management, and physical activity/sedentary behavior. In one of the first articles on this topic, Puterman et al. (2010) reported that exercise might have a great effect on attenuating telomere attrition. In a larger study, Sjogren et al. (2014) showed that avoiding sedentary behavior might also attenuate telomere attrition.
It seems that avoiding sedentary behavior leads to well-being, but adding exercise daily may lead to even better well-being in terms of telomere preservations. Elite athletes have longer TL than their counterparts (Muniesa et al., 2017) and master athletes (around 40-50 years old) also have longer TL than age-matched non-athletes (Østhus et al., 2012;Simoes et al., 2017).
The mechanisms by which exercise could attenuate cell senescence and preserve TL have being studied especially in animal models (Werner et al., 2009;Cunha et al., 2018). Werner et al. (2009) reported that increased levels of PA are associated with elevated telomerase activity and suppression of several apoptosis proteins, such as p53 and p16. On the other hand, Cunha et al. (2018) showed that more vigorous exercise (above lactate threshold) could also diminish p53 activity, but also decrease the expression of proteins related to telomere protection (sheltering proteins). These authors suggest that this may happen because in a more balanced environment, there is no need to a larger expression of sheltering proteins.
In summary, the body of available evidence shows that one of the factors that exercise may lead to an attenuated telomere attrition and aging is the maintenance of good body composition (Sousa-Victor et al., 2015;Simoes et al., 2017). Which, in turn, it likely helps to maintain a good metabolic balance, as well as healthy states of oxidative stress and inflammatory status (von Zglinicki, 2002).

THE INTERPLAY
Oxidative stress seems to play a central role in which TL could increase the cancer risk. The Long Island Breast Cancer Study Project (Gammon et al., 2002) explored the relationship between TL and breast cancer risk, examining whether dietary intake of antioxidants would have any effect. The authors of this project reported that a moderate increase in breast cancer risk was observed among women with the shortest telomeres and with lower dietary and supplemental intake of beta-carotene, vitamin C or E. These results provided the evidence that breast cancer risk may be affected by TL among premenopausal women or women with low dietary intake of antioxidants or in an unbalanced oxidative stress status (Shen et al., 2010). Antioxidants dietary intake (through supplementation or general food) is important to a healthier redox status and TL maintenance (Noorimofrad and Ebrahim, 2018). However, the regulation of endogenous antioxidant capacity by regular exercise has led to favorable results (Radak et al., 2008;Sousa et al., 2017).
In that context, exercise seems to play a central role in cancer treatment and prevention, and also in TL maintenance. Nevertheless, the available literature shows that this may occur by different pathways. Biologically, longer TL means that a cell is young and has not been divided too many times yet. Roughly, cell division happens when another cell is needed (Allsopp et al., 1995). If the need for new cells is greater than the capacity of cell to multiply themselves, the tissue starts to lose their function little by little, that is aging (Allsopp et al., 1995). However, cancer cell function as the opposite of that, this mutant genes have a high division rate and normally have longer TL and telomerase activity (Macieira-Coelho, 1994;Brinkley, 2001).
Rapidly dividing cells causes cancer, and some of those (i.e., esophageal and ovarian) are associated with shorten TL (Risques et al., 2007;Wentzensen et al., 2011). McGrath et al. (2007 observed, in both men and women, that those with bladder cancer had shorter TL than controls. They also evaluated the effect of smoking on TL and found significantly shorter telomeres in healthy individuals who smoked than in those who did not (McGrath et al., 2007). A series of epidemiological studies found that having shorter overall TL from blood leucocytes was significantly associated with an increased risk of bladder, head and neck, lung, and renal cell cancers (Zheng et al., 2010). Pooley et al. (2010) found a strong association between shorter telomere length and breast cancer status in the retrospective study, with an OR of ∼15 for women in the bottom quartile of TL compared to the top quartile (Pooley et al., 2010). They also stated that the prospective breast cancer study showed evidence of an association, but the effect was not statistically significant and the majority of telomere attrition occurred after cancer diagnosis rather than before or during cancer (Pooley  Martinez-Delgado et al. (2012) reported that shorter telomere length is associated with increased ovarian cancer risk in familial and sporadic cases (Martinez-Delgado et al., 2012). Patients with colorectal cancer (CRC) displayed clear evidence of telomere attrition compared with controls (Maxwell et al., 2011). Furthermore telomere length could be a non-invasive blood biomarker to pre-screen for risk of advanced adenomas (AAs); PBL telomeres are shortened in patients with colorectal neoplasia (Riegert-Johnson et al., 2012). Short TL and subsequent genomic instability contribute to malignant transformation, and it is therefore contingent that people with minor TL are at higher risk for different types of cancer (Martinez-Delgado et al., 2012).
Furthermore, Ma et al. (2011) performed a meta-analysis and reported that, although shortened TL are associated with increase risk of bladder and lung cancer, this association was not found for breast cancer. On the other hand, another study showed that an increased breast cancer mortality following diagnosis to women with elevated telomere shortening (Duggan et al., 2014). More recently, Sanft et al. (2018) performed a randomized controlled trial in women with breast cancer, from the total sample size (n = 151), 93 were submitted to a lifestyle intervention that included weight loss strategies with exercises and nutrition and 58 to usual care. The authors reported that, the lifestyle intervention in breast cancer may lead to TL maintenance, but the effect may be limited to patients in the earlier stages (0/I) of breast cancer.
Taken separately, in a meta-analysis was showed evidence that perceived psychological stress is associated with a decrease in TL (Mathur et al., 2016). In several other observational studies, was evidenced that exercise as could reduce telomere attrition and contribute to their length maintenance (Kim et al., 2012;Ornish et al., 2013;Puterman et al., 2018). And also that most types of cancer, although are closely related to telomere, they are also strongly associated with lifestyle (Eisenberg, 2011;Kerr et al., 2017). Therefore, it is reasonable to suggest that shorter telomere length may not be a direct cause for cancer, but a result of inadequate lifestyle choices that possibly led to both increased risk of cancer and shorter TL (Figure 1). The possible pathway between cancer and TL still needs clarification.

CONCLUSION
In summary, there is evidence that exercise leads to less telomere attrition and exercise also may diminish the risk of cancer, these two outcomes are possible intermediated by a reduction in oxidative stress and chronic inflammation. Although, there is evidence that shorter TL are associated with cancer, the possible mechanisms that one may lead to the other remains to be clarified. We assume that humans under cancer treatment may suffer a great decrease in quality of life, which may increase sedentary behavior and lead to increased telomere attrition. And those humans with already shorter TL likely lived under a poor lifestyle and might have an increased risk to have cancer.

AUTHOR CONTRIBUTIONS
NN designed the review and collected the data. All authors contributed by writing, editing, and reviewing the manuscript, and read and approved the final manuscript.

ACKNOWLEDGMENTS
Special thanks to Dionysios Chorafas MSc., for the linguistic editing of the manuscript.


extfieldsofstudy: Biology Medicine

text:
Complex rearrangements and oncogene amplifications revealed by long-read DNA and RNA sequencing of a breast cancer cell line

The SK-BR-3 cell line is one of the most important models for HER2+ breast cancers, which affect one in five breast cancer patients. SK-BR-3 is known to be highly rearranged, although much of the variation is in complex and repetitive regions that may be underreported. Addressing this, we sequenced SK-BR-3 using long-read single molecule sequencing from Pacific Biosciences and develop one of the most detailed maps of structural variations (SVs) in a cancer genome available, with nearly 20,000 variants present, most of which were missed by short-read sequencing. Surrounding the important ERBB2 oncogene (also known as HER2), we discover a complex sequence of nested duplications and translocations, suggesting a punctuated progression. Full-length transcriptome sequencing further revealed several novel gene fusions within the nested genomic variants. Combining long-read genome and transcriptome sequencing enables an in-depth analysis of how SVs disrupt the genome and sheds new light on the complex mechanisms involved in cancer genome evolution.

Genomic instability is one of the hallmarks of cancer, leading to widespread copy number variations, chromosomal fusions, and other sequence variations (Hanahan and Weinberg 2011). Structural variations, including insertions, deletions, duplications, inversions, or translocations at least 50 bp in size, are especially important to cancer development, as they can create gene fusions, amplify oncogenes, delete tumor suppressor genes, or cause other critical changes to contribute to the evolution of a cancer genome (Mitelman et al. 2007). Detecting and interpreting these structural variations is therefore a crucial challenge as we try to gain a more complete understanding of cancer genomes (Nik-Zainal et al. 2012).
Cancer genomics has been greatly aided by the advances in DNA sequencing technologies over the last 10 yr (Watson et al. 2013). The first whole-genome analysis of a cancer genome was reported in 2008 (Ley et al. 2008), and today large-scale efforts such as The Cancer Genome Atlas (Kandoth et al. 2013) or the International Cancer Genome Consortium (The International Cancer Genome Consortium 2010) have sequenced thousands of samples using short-read sequencing to detect and analyze commonly occurring mutations, especially single nucleotide and other small variations. However, these projects have performed somewhat limited analysis of structural variations, as both the false positive rate and the false negative rate for detecting structural variants from short reads are reported to be 50% or more Huddleston et al. 2017). While short-read sequencing has no doubt revolutionized cancer genomics, this latter observation is troubling. Furthermore, the variations that are detected are rarely close enough to determine whether they occur in phase on the same molecule, limiting the analysis of how the overall chromosome structure has been altered.
Addressing this critical void, we sequenced the HER2 (ERBB2)amplified breast cancer cell line SK-BR-3 using long-read sequencing from Pacific Biosciences. SK-BR-3 is one of the most widely studied breast cancer cell lines, with applications ranging from basic to preclinical research (Lewis Phillips et al. 2008;Navin et al. 2011;Ichikawa et al. 2012). SK-BR-3 was chosen for this study due to its importance as a basic research model for cancer and because the SK-BR-3 genome contains many of the common features of cancer alterations including a number of gene fusions, oncogene amplifications, and extensive rearrangements. Critically, the amplifications and genome complexity observed in SK-  have been demonstrated to be representative of patient tissues as well (Neve et al. 2006).
Taking full advantage of the benefits of the new long-read sequencing technology, we applied a split-read and within-read mapping approach to detect variants of different types and sizes. This allows us to develop a comprehensive map of structural variations in the cancer and study for the first time how the rearrangements have occurred with base-pair level accuracy. Furthermore, combining genomic variant discovery with Iso-Seq full-length transcriptome sequencing, we discover new isoforms and characterize several novel gene fusions, including some that required the fusion of three separate chromosome regions. Finally, using the reliable mapping and coverage information from long-read sequencing, we show that we can reconstruct the progression of rearrangements resulting in the amplification of the ERBB2 oncogene, including a previously unrecognized inverted duplication spanning a large portion of the region. Using long-read sequencing, we document a great variety of mutations including complex variants and gene fusions far beyond what is possible with alternative approaches.

Results
We sequenced the genome of SK-BR-3 using Pacific Biosciences (PacBio) SMRT long-read sequencing (Eid et al. 2009) to 71.9× coverage (based on the reference genome size) with an average readlength of 9.8 kb (Supplemental Fig. S1). For comparison, we also sequenced the genome using short-read Illumina paired-end and mate-pair sequencing to similar amounts of coverage. To investigate the relevant performance of long and short reads for cancer genome analysis, we perform an array of comparisons in parallel using both technologies.

Read mapping and copy number analysis
Long reads have more information to uniquely align to the genome than short reads do, resulting in overall better mapping qualities for long reads (Supplemental Fig. S2; Lee and Schatz 2012). Using BWA-MEM (Li 2013) to align both data sets, 69% of Illumina short paired-end reads (101-bp reads, 550-bp fragment length) align with a mapping quality of 60 compared to 91.61% of reads from the PacBio long-read sequencing library (Supplemental Fig. S2; Supplemental Table S1). We also observed a smaller GC bias in the PacBio sequencing compared to the Illumina sequence data which enables more robust copy number analysis and generally better variant detection overall (Supplemental Fig. S3).
The average aligned read depth of the PacBio data set across the genome is 54×, although there is a broad variance in coverage attributed to the highly aneuploid nature of the cell line (Supplemental Fig. S4). The short reads showed a few regions of extreme amplification (>100-fold) that were not detected by the long reads, although subsequent analysis showed these regions were highly enriched for low mappability regions (Dolgalev et al. 2017) in the genome and therefore most likely to be mapping artifacts (Supplemental Fig. S5). Using the long-read alignments, we segmented the genome into 4083 segments of different copy number states with an average segment length of 747.0 kbp. The unamplified chromosomal regions show an average coverage of 28×, which we consider the diploid baseline for this analysis. Thus, the average copy number is approximately twice the diploid level, which is consistent with previous results characterizing SK-BR-3 as tetra-ploid on average (Navin et al. 2011), and with any given locus being heterogeneous in copy number across the cell population.
Assuming a diploid baseline of 28×, the locus spanning the important ERBB2 oncogene (17q12) is one of the most amplified regions of the genome with an average of 33.6 copies (average read coverage of 470×). A few other regions show even greater copy number amplification, including the region surrounding MYC, with 38 copies. Other oncogenes are also amplified, with EGFR at seven copies and BCAS1 at 16.8 copies, while TPD52 lies in the middle of an amplification hotspot on Chromosome (Chr) 8 and is spread across eight segments with an average copy number of 24.8. The locus 8q24.12 containing the SNTB1 gene is the most amplified region of the genome with 69.2 copies (969× read coverage). In addition to being the most amplified protein-coding gene in this cell line, SNTB1 is also involved in a complex gene fusion with the KLHDC2 gene on Chr 14 (see below). Copy number amplifications are distributed throughout the genome across all chromosomes (Supplemental Fig. S5). Every chromosome has at least one segment that is tetraploid or higher, and these amplified regions account for about one third (1.07 Gbp) of the genome. Extreme copy number amplifications, above 10-ploid (>140× coverage), appear on 15 different chromosomes for a total of 61.1 Mbp, with half on Chr 8 (30.1 Mbp). There is a total of 21.3 Mbp of 20-ploid sequences across five chromosomes, with 20.0 Mbp on Chr 8 and 1.3 Mbp distributed across Chromosomes 17, 7, 21, and 1. In addition to containing the greatest number of base pairs of 20-ploid sequence, Chr 8 also has 101 segments of 20-ploid sequence compared to only four total segments from Chromosomes 7, 16, 17, and 21. Chr 8 thus has far higher levels of extreme copy number amplification than all other chromosomes combined.

Structural variant analysis
We aligned the long reads to the reference using NGMLR (Sedlazeck et al. 2018), a read mapping algorithm optimized for long single-molecule sequencing reads, and analyzed the alignments for structural variations using Sniffles (Sedlazeck et al. 2018). Sniffles was specifically designed for long-read SV analysis and identifies them from both split-read and within-read alignments requiring at least 10 split reads to call a structural variant. Sniffles found a total of 76,776 variants that were 10 bp or larger, and of these, 17,313 variants were structural variants (50 bp or larger), composed of 8909 (51%) insertions, 6947 (40%) deletions, 1018 (6%) duplications, 279 (2%) inversions, and <1% total of translocations and special combined variant types (Fig. 1). Our work with several other genomes shows that the vast majority of variants called using this combination of algorithms are correct (also see below) (Sedlazeck et al. 2018). Within all of the variants detected using the long reads, 1725 variants intersect transcribed regions of 361 of the 616 genes in the COSMIC Cancer Gene Census (Futreal et al. 2004), and 172 of these genes are hit by structural variants with a minimum size of 50 bp (Supplemental Tables  S4, S8). Counting only sequences identified in GENCODE as exons, a total of 58 variants intersect the exons of 46 different Cancer Gene Census genes, including breast cancer genes such as APOBEC3B and CDH1. The deletion in APOBEC3B is consistent with a germline variant previously observed to fuse APOBEC3A and APOBEC3B, which is suggested to confer increased cancer susceptibility (Nik-Zainal et al. 2014).
For comparison, we also called structural variants from a standard paired-end short-read sequencing library, using our SURVIVOR algorithm (Jeffares et al. 2017) to form a high-quality consensus call set from three different short-read variant callers (Manta [Chen et al. 2016], DELLY [Rausch et al. 2012], and LUMPY [Layer et al. 2014]), requiring that at least two of these variant callers identified the same variant. We have found this approach reduces the false positive rate without substantially reducing sensitivity (Jeffares et al. 2017). The total number of short-read structural variants in the consensus set was 4174, composed of 2481 (59%) deletions, 603 (14%) translocations, 580 (14%) inversions, 448 (10%) duplications, and 62 insertions (1.4%). Comparing the counts, the short-read consensus has a much smaller number (only 24%) of total variants than the long-read set (see Supplemental Figs. S9-S16 for examples of variants not detected by the short-read analysis). This difference is largely driven by the lack of insertions in the short-read call sets: DELLY and Manta report a small number of insertions, but LUMPY does not attempt to report any (Rausch et al. 2012;Layer et al. 2014;Chen et al. 2016). To further address the limited number of insertions detected, we also ran a new insertion-finding algorithm called PopIns (Kehr et al. 2016). Using the recommended settings, PopIns finds 579 insertions that it could anchor to the genome, of which six were also found by DELLY and 11 were also found by Manta (two were found by both). Overall, this raised the total number of insertions detected by two or more shortread mapping algorithms to 77 (Figs. 1B, 2B). Supplemental Table S5 shows the count of variant calls of each type, and Figure 1C shows the counts for variants that are at least 1 kbp in size. For noninsertions, we also note that the short-read SV callers are highly enriched for false positive calls, especially false translocations (see below). The disagreement between the short-and long-read variant calls does not appear to be related to coverage. Using SAMtools (Li et al. 2009) bedcov, we found the mean short-read coverage for variants found by both short and long reads was 25.6× coverage at the breakpoint, while the short-read coverage for variants only detected by long reads was 51×.
In parallel, we performed assembly-based variant calling with Assemblytics (Nattestad and Schatz 2016) using a de novo assembly of the long reads using FALCON (Chin et al. 2016) and the short reads using ALLPATHS-LG (Supplemental Note S1; Gnerre et al. 2011). The long-read assembly achieved a 2.4-Mbp contig N50 and showed good sensitivity for many structural variant types, especially insertions and deletions <1 kbp in size (Supplemental Figs. S6-S8) and found variations within many hundreds of ALU sequences (Supplemental Table S3). However, we also found this approach misses many long-range variants due to splitting of the assembly graph into contigs at or near the branch points caused by large variants. Therefore, the long-range variants are not always captured well within the contigs, and the evidence can be skewed to the end of the contigs where alignment and assembler errors are more common. Consequently, assemblybased variant calling is therefore not ideal for this class of variants. The short-read assembly was even more limited, as the contig N50 was only 3.2 kbp, and only a small number of SVs could be found (Supplemental Tables S2, S3).
Our initial expectation was that approximately the same number of insertions and deletions would be present due to normal human genetic variation. However, the long-read variant call set has a ratio of 1.28:1 insertions to deletions. This insertional bias has been seen previously and suggests an underestimate of the lengths of low-complexity regions in the human reference genome (Chaisson et al. 2015). In support of this analysis, using the repeat annotation tracks from the UCSC Genome Browser, we found that 52 Sniffles insertions are within annotated microsatellites, 5015 insertions are within simple repeats, and 6027 are within regions identified by RepeatMasker (Smit et al. 2013(Smit et al. -2015, in agreement with the prior studies. As long-range variants are of particular interest in cancer genomics, we performed several analyses specific to this subset of variants. We define long-range variants as those that are either (1) between different chromosomes, (2) connecting breakpoints at least 10 kbp apart within the same chromosome, or (3) inverted duplications. These long-range variants indicate novel adjacencies joining chromosomal regions that were originally distant in the genome. This causes novel sequences to be formed at the junction, potentially leading to gene fusions, large deletions or duplications, and other aberrant genomic features. Split reads provide a robust signal for detecting these long-range variants and chromosomal rearrangements. Within the long-read Sniffles call set, we found 665 variants in this long-range class ( Fig. 1A; Supplemental Table  S6), 125 of which were between different chromosomes. From the SURVIVOR short-read consensus calls, 1493 are long-range variants, with 603 of these being between different chromosomes.
Focusing on the long-range variants, we analyzed the intersections between the Sniffles and SURVIVOR (two-caller consensus) call sets. Compared to the SURVIVOR consensus call set, Sniffles detects the same 461 and an additional 204 variants, whereas the short-read SURVIVOR consensus detects an additional 1032 ( Fig. 2A). We randomly selected 100 variants from each sub-set for PCR plus Sanger validation, with 100 calls from the intersect, 100 Sniffles calls not shared by SURVIVOR (two-caller consensus), and 100 SURVIVOR 2-calls not shared by Sniffles. Within each randomly selected group of 100 variants, some variant calls could not be validated due to primer issues or other technical issues, so the final validation rates are calculated as successful Sanger validation counts out of the total valid attempts. As expected, the variants called by both Sniffles and SURVIVOR had the highest validation rate of 82.8% (77/93) (Supplemental Table  S7). Of the 16 variants called by both short-and long-read approaches that failed validation, five calls were translocations, and the rest were paracentric. All of the reported variants had at least 10 supporting long reads (the minimum threshold we used for Sniffles), and ranged from 10 to 31 supporting reads (mean: 17.6×), except for one outlier with 113 supporting reads. Seven of the PCR attempts produced no recognizable product even after multiple attempts and two primer designs. The other nine attempts produced a weak product or a multibanded product but failed to report the expected sequence during Sanger validation, and so we consider them as failing validation. Given the strong support from both PacBio and Illumina sequencing, we attribute the failures as either inadequate primer design and/or other systematic errors in the validation protocol. Of the calls unique to one method, long-read Sniffles variants have a validation rate more than twice that of the short-read variants: 48.2% (26/54) compared to 21.3% (17/80). Furthermore, extrapolating the validation rates for these subsets, the overall validation rate of Sniffles calls is 72%, while the SURVIVOR 2-calls is only 29.6%. We emphasize that this is the validation rate for the most complicated long-range variants present in the genome, and our work with structural variant detection in other long-read data sets reached 94% to over 99% (Sedlazeck et al. 2018). We also note that several of these long-range variants were previously found through RNA-seq and confirmed with mate-pairs in previous studies (Table 1; Kim and Salzberg 2011).
Further supporting this higher validation rate of long-read variants, the Sniffles variants were also more likely to occur at the breakpoint of a copy number variant than their short-read counterparts. Specifically, 58.3% of the Sniffles unique variants show a matching copy number variant, compared to only 23.2% of the SURVIVOR unique consensus variants, where 58.1% of the variants shared by both sets show a matching CNV (Fig. 2B,C). Similar results were also found using the short reads for segmentation. The high rate of CNV matching for the shared set indicates that copy number evidence can serve as a measure of confidence in a variant call. CNV matching provides additional support for the majority of the Sniffles unique calls, though it does not exclude others that may be copy number neutral variants such as balanced translocations. The low rate of CNV support for the short-read consensus suggests that a larger proportion of these variants are either false positives or Sniffles is not sensitive enough to capture them. Reducing the threshold in Sniffles to five split reads (instead of the 10 split reads employed throughout this analysis) captures another 134 of the short-read consensus variants out of 1032, so there appears to be little long-read evidence of these additional variants.
Characterization of the ERBB2 copy number amplification 8 are clustered in three major hotspots. The ERBB2 oncogene, originally located on Chr 17, is amplified to, on average, 32.8 copies, while most of the remainder of Chr 17 is present in just two copies, consistent with selection against gains of tumor suppressor proteins on Chr 17, such as TP53 and BRCA1. The amplified region that includes ERBB2 contains five translocations (Fig. 3A) into the hotspot regions on Chr 8, as well as an inverted duplication. Each of the six variants mark the site of an abrupt change in copy number, meaning they mediated the overall amplification to 32.8-fold, and all six were validated by directed PCR and Sanger sequencing. It is notable that the inverted duplication was not identified by any of the short-read variant callers, although it is clearly visible in the long-read alignments and is automatically identified by Sniffles.
The ERBB2 oncogene appears to have been amplified to such a great extent due to its association with the highly mutated hotspots in Chr 8 and suggests a remarkably complex and punctuated mutational history (Fig. 3B). The long-range variants within the amplified region containing ERBB2 were studied to determine whether the number of split and reference-spanning reads at each breakpoint are consistent with the copy number profile, which was found to be true for all five translocations and the inverted duplication. In order to determine the order in which these six events took place, we derived the most parsimonious reconstruction factoring in a couple of important assumptions established within population genetics, analogous to the widely used infinite sites model used in population genetics (Kimura 1969). First, we assume that variants we observe have taken place only once, since it is extremely unlikely that the same long-range variant at the same two breakpoints would recur down to base-pair resolution. Second, once a variant has occurred and created an observable breakpoint, the breakpoint would not be repaired in some copies of the sequence and not others. Therefore, all reference-spanning reads represent an ancestral state and not a repaired breakpoint. In this analysis, the long-range variants have more reliability to reconstruct the genomic history rather than SNPs because those two simplifying assumptions are extremely unlikely to be violated when two breakpoints are involved.
Given these conditions, we conclude that the orange segment (A-F) in Figure 3B must have translocated first, as the other breakpoints are shared on the leftmost edge (variant A). Next, the yellow segment shown in Figure 3B derived from the orange segment because otherwise, variant A must have occurred more than once. Applying the same logic, the green segment must have derived from the yellow segment because it shares variant E, and it is not yellow derived from green because that would violate assumption 2 by requiring that variants C and D were repaired. Variants C and D appear to co-occur in the same sequences because the copy number is the same between those two parts of the green segment and because the other sides of the variants are at breakpoints within only 1.5 Mb of each other. The only uncertainty in the ordering of events is that the purple segment could have derived from any of the segments sharing variant A: the orange, yellow, or green segments. There is not enough information to determine which of these segments it came out of, but we can conclude that it only came out of one of them, given assumption 1 that precludes multiple occurrences of the same variant.

Complex gene fusions captured fully by long reads
In addition to genome sequencing, we performed long-read transcriptome sequencing using PacBio Iso-Seq to capture full-length transcripts. Although traditional short-read RNA-seq approaches allow isoform quantification, in many cases these reads are too short to reconstruct all isoforms, even with paired-end analysis, exon abundance, or other indirect measurements. Instead, long reads overcome such limitations by spanning multiple exon junctions and often covering complete transcripts. This makes it possible to exactly resolve complex isoforms and identify large transcripts, without the need for statistical inference (Sharon et al. 2013;Weirather et al. 2015;Wang et al. 2016). Iso-Seq reads were consolidated into isoforms using the SMRTAnalysis Iso-Seq pipeline (Gordon et al. 2015). In total, 1,692,379 isoforms (95.7%) mapped uniquely to the reference genome. The Iso-Seq RNA sequence reads indicated a total of 53 putative gene fusions each with at least five Iso-Seq reads of evidence (Supplemental Table S9). We further refined this candidate set using SplitThreader (Nattestad et al. 2016a) to exclude variants not supported by genomic structural variations, especially to account for any residual sequencing error or mapping errors in the data. Specifically, SplitThreader searches for a path of structural variations linking the pair of genes in the putative gene fusion, requiring that the variants bring the genes together within a 1-Mbp distance. Out of 53 candidate gene fusions, SplitThreader found genomic evidence for 39 of these: 15 are the high-quality gene fusions with a genomic path between the gene bodies of, at most, 10 kbp, shown in Table 1; 19 fusions overlap the first 15 (sharing the same variant and often one of the genes), and five fusions (three nonoverlapping) have paths longer than 10 kbp, leaving 14 candidate gene fusions with no genomic paths.
Three of the gene fusions had no single variant directly linking the genes, but SplitThreader discovered that the genes could be linked by a series of two or even three variants. One of these, CPNE1-PREX1 had been discovered previously using RNA-seq data and validated using genomic PCR as a two-variant gene fusion (Chen et al. 2013). We have now confirmed this by showing long reads that not only capture the two variants but capture them together in a single read along with robust alignments to both genes (Supplemental Fig. S19). CYTH1-EIF3H had been discovered previ-ously with RNA-seq and been validated with RT-PCR (Edgren et al. 2011), but it was not known to be a "2-hop" gene fusion (taking place through a series of two variants) until now. This fusion was also captured in full by several individual SMRT-seq reads that contain both variants and have alignments in both genes (Supplemental Fig. S18). Interestingly, we discovered a novel 3-hop gene fusion between KLHDC2 and SNTB1, which has been misreported before as only taking place through two variants (Asmann et al. 2011). We observe both the previously reported 2-hop path (600,326 bp) and this additional 3-hop path (9837 bp), which would both result in the same gene fusion. Given the shorter distance for the 3-hop gene fusion, we were able to find direct linking evidence for the 3-hop fusion between these two genes. We observe 37 reads that stretch from one gene to the other through all three variants, bringing the genes within a distance of just 9837 bp across three different chromosomes ( Fig. 4; Supplemental  Fig. S17). Due to the long distance between the genes through the previously reported 2-hop fusion, we believe the 3-hop fusion is more likely to produce the observed fusion transcript.
Most of the gene fusions observed are contained within a few of the most rearranged chromosomes. Four gene fusions take place within Chr 20, which is rich in intra-chromosomal variants, while Chr 8 is involved in six gene fusions both intra-and inter-chromosomally. The genomic variant fusing TATDN1 and GSDMB is one of the variants contributing to the amplification of the ERBB2 oncogene. All of the gene fusions are captured fully with individual SMRT-seq reads that align to both genes, with some novel variants affecting important cancer genes (e.g., PVT1 and RAD51B). See long-read alignments spanning all 15 gene fusions in Supplemental Note S2 and Supplemental Figures S20-S32.

Discussion
Advances in long-read sequencing have produced a resurgence of reference quality genome assemblies and exposed previously B A Figure 3. Reconstruction of the copy number amplification of the ERBB2 oncogene. (A) Copy number and translocations for the amplified region on Chr 17 that includes ERBB2 showing the relations to Chr 8. Note Chr 8 has extensive rearrangements shown by the green intra-chromosomal arcs. (B) Sequence of events that best explains the copy number and translocations found in this region. Segment 1 (orange) first translocated into Chr 8, followed by the segment 2 (yellow) translocating to a different place on Chr 8. Then, the segment 3 (green) was duplicated from segment 2 by an inversion of the piece between variants D and E along with a 1.5-Mb piece of Chr 8 that was attached at variant E, all of which then attached at variant C. The whole green segment including the 1.5 Mb of Chr 8 then underwent an inverted duplication at variant D. The purple segment could have come from the orange, yellow, or green sequences since it only shares breakpoint A. Additionally, there is a deletion of 10,305 bp between breakpoints D and E.
hidden genomic variation in healthy human genomes (Chaisson et al. 2015;Pendleton et al. 2015;Seo et al. 2016). Now, we have applied long-read sequencing to explore the hidden variation in a cancer genome and have discovered nearly 20,000 structural variations present, most of which cannot be found using short-read sequencing and many are intersecting known cancer genes. More than twice as many of the copy number amplifications could be explained through long-range variants identified by long-read sequencing compared to short-read sequencing. We further found the ERBB2 oncogene to be amplified through a complex series of events initiated by a large translocation into the highly rearranged hotspots of Chr 8, where the sequence was then copied dozens of times more, with further translocations and inverted duplications resolved only by the long reads. Furthermore, we find 20 additional inverted duplications throughout the genome, highlighting the importance of this underreported structural variation type. Overall, using long-read sequencing, we see that far more bases in the genome are affected by structural variation compared to SNPs.
Using long-read transcriptome sequencing, we capture full gene fusion isoforms, and by combining this with our genomic variant discovery, we discover several novel gene fusions in this seemingly well-characterized cell line. Notably, we uncover for the first time a gene fusion that takes place through a series of three variants: KLHDC2-SNTB1 through the fusions of Chromosomes 8, 14, and 17, captured fully by 37 genomic SMRT-seq reads. In a single cancer genome, we discovered three gene fusions that take place through a series of two or more variants, suggesting that such multi-hop gene fusions could also be common in other cancers, although they will be exceedingly difficult to discover using short-read sequencing. Conducting a similar search for multi-hop gene fusions in other highly rearranged cancers could reveal other instances of a complex type of variation.
The differences between variants found with long reads versus short reads is likely due to an interplay of sequencing technology and algorithmic approach. The primary advantages of long reads are better mapping through repetitive elements that often flank SVs and an increase in the probability that a SV breakpoint will be spanned by individual reads. These advantages are offset by the increased error rate that makes them more difficult to map and analyze, although new mapping and SV detection tools are now available that can largely overcome these challenges. For short-read analysis, using paired-end or mate-pair sequencing can partially offset the short-read lengths, but they still have relatively poor sensitivity using current approaches. In addition to the results presented here, a previous study by Hillmer et al. (2011) analyzed SVs in SK-BR-3 using long-range mate-pairs. In this work, they generated 68.4× physical coverage of mate-pairs averaging 8.2 kbp and yet only found 1145 SVs, most of which were deletions (606), followed by inversions (191), and other complex intra-chromosomal rearrangements (158). This represents <10% of the variants that we could detect using long reads. Further algorithmic advances may be possible to improve accuracy, and we currently recommend using a consensus approach for short-read analysis to alleviate false positives. It may also be possible to improve sensitivity of certain types of variants using focused methods. For example, ARC-SV (Arthur et al. 2018) and SVelter (Zhao et al. 2016) were recently developed to focus on distal and inverted duplications from short-read sequencing, although these classes of variations are a minority in this sample. What is most needed is a robust method to detect insertions from short reads, as they are currently A B Figure 4. The KLHDC2-SNTB1 gene fusion in SK-BR-3 occurs through a series of three variants and is directly observed to link the two genes in several individual SMRT-seq reads (A), one of which is shown in detail in B. not well captured by either mapping-based or assembly-based approaches.
We have shown that long-read sequencing can expose complex variants with great certainty and context, suggesting that more multi-hop gene fusions, inverted duplications, and complex events may be found in other cancer genomes. Having observed complex variants such as inverted duplications with the increased informational context of long reads, the resulting variant signatures could make these events more observable even using standard short-read sequencing. However, there may be many other types of complex variations present in other cancer genomes that were not found in SK-BR-3, so it is essential to continue building a catalog of these variant types using the best available technologies. Long-read sequencing is an invaluable resource to capture the complexity of structural variations on both the genomic and transcriptomic levels, and we anticipate widespread adoption for research and clinical practice as the costs further decline.

Methods Sequencing
Long-read sequencing was performed using the Pacific Biosciences Single-Molecule Real-Time (SMRT) sequencing technology with P6C4 SMRT cell chemistry. After selecting the longest subread from each polymerase PacBio read, our sequencing of SK-BR-3 yielded a mean read length of 9872 bp, where the longest read was 71,518 bp. Total coverage of the genome is 71.9× (79.0× if redundant sequences from the same polymerase reads are included) where × refers to the number of reads that cover the average genomic base. The coverage of reads at least 10 kbp long is 51.0×, and the coverage of reads at least 20 kbp long is 13.3×. These read depth values and those in Supplemental Figure S1B are based on a female genome size of 3,101,804,739 bp, the total lengths of Chromosomes 1-22 and Chr X in hg19.
For short-read variant calling, Illumina sequencing was performed on a 550-bp paired-end library (2 × 250 bp). This library produced a total of 795,942,102 reads and 64.2× genome coverage based on the same female genome size. For the short-read assembly, Illumina sequencing was performed on 180-bp paired-end overlapping library (2 × 100-bp reads), as well as 2-3 kbp and 5-10 kbp mate-pair libraries.

Alignment and variant calling
The hg19 reference genome (the 1000 Genomes version) was used for all analysis. Results aligning to GRCh38 are expected to be similar as no major differences were introduced in the ERBB2 locus and only a minority of the bases changed genome-wide. Reads were aligned to the reference using NGMLR (v0.2.1) (Sedlazeck et al. 2018), and Sniffles (v1.0.6) (Sedlazeck et al. 2018) was used to call variants from long-read split alignments using the recommended parameters. Variants were called on the short-read variant-calling Illumina sequencing data set using Manta (Chen et al. 2016), DELLY (Rausch et al. 2012), LUMPY (Layer et al. 2014), and PopIns (Kehr et al. 2016), and a consensus was taken using SURVIVOR (Jeffares et al. 2017) with the recommended parameters, requiring two of these variant callers to support the same variant , except where noted otherwise. Because there is often variability in the reported location of SVs, this process allows for merging of SVs that have breakpoints within 1 kbp of each other, as long as the type of variant is the same. For purposes of intersecting Sniffles and SURVIVOR long-range variant call sets, for which breakpoints must be 10 kbp apart, we used BEDTools (Quinlan and Hall 2010) pair-to-pair with a slop parameter of 1000 to test if variants shared both breakpoints with matching strands within a 1000-bp range. Copy number segmentation was computed using SplitThreader (Nattestad et al. 2016a), which internally uses the DNAcopy R package for circular binary segmentation (Olshen et al. 2004). The SplitThreader source code is also available as a Supplemental File. The Circos plot in Figure 1A was generated using Circa (http://omgenomics.com/circa). Cancer gene intersects were determined using BEDTools pairtobed to intersect Sniffles variants down to 10 bp in size with the GENCODE hg19 annotation (Harrow et al. 2012) and filtered by matches to the COSMIC Cancer Gene Census (Futreal et al. 2004).

Mapping comparison
In order to compare the mappability of long and short reads, we aligned both the paired-end Illumina sequencing and the PacBio long-read sequencing data sets to the hg19 reference genome using BWA-MEM (Li 2013). The Illumina sequencing was performed using a 550-bp paired-end library with each read being ∼250 bp of sequence. We trimmed these reads to 101 bp and compared both of these against the PacBio data set. All three read sets were aligned using default parameters, except that the PacBio reads were aligned using the pacbio alignment mode in BWA-MEM (-x pacbio). The maximum mapping quality in BWA-MEM is 60, and the minimum is 0. Using the same aligner allows us to better compare mapping quality scores for the reads. We analyzed the mapping quality from each type of sequencing in two different ways, by individual reads and by binned windows in the genome. First, we selected the best alignment by mapping quality for each read and counted the number of reads in each category: mapping quality of 60, mapping quality between 1 and 59, mapping quality of 0, or unmapped. Alignment of PacBio sequence reads resulted in 91.6% of reads mapping with a mapping quality of 60, compared to only 71.2% of Illumina reads (69.0% of the 101-bp trimmed reads). A greater fraction of reads from PacBio long-read sequencing map uniquely to the genome compared to short reads from Illumina sequencing (Supplemental Fig. S2; Supplemental Table S1).
In order to determine the effect of GC content (the fraction of guanine and cytosine as opposed to adenine and thymine in a particular region), we counted the GC-fraction of each 10-kbp window in the genome, excluding those containing Ns in the reference, and calculated the read coverage from each data set. The read depth of each 10-kb bin is shown in Supplementary Figure S3 on a log scale versus the GC fraction, along with a Lowess fit for each data set. There is a higher GC-bias in the Illumina data sets compared to the PacBio data set, as seen by a lower read depth in bins with a higher GC fraction, while for SMRT sequencing there is a much lower bias.
To determine the read depth per chromosome, we used BEDTools to find the distribution of read depth for each chromosome for the PacBio, Illumina 250-bp, and Illumina 101-bp data sets. These are shown as a violin plot of Gaussian kernel distributions for each chromosome in Supplemental Figure S4. The shapes of the distributions are largely consistent between sequencing technologies.

Iso-Seq and gene fusion analysis
PacBio Iso-Seq sequencing was performed in four size batches (0.8-2, 2-3, 3-5, and 5-10 kb). The Iso-Seq data were processed using the SMRTAnalysis (version 2.3) Iso-Seq pipeline, which generated 441,932 high-quality (HQ), full-length Quivered consensus sequences, which were then aligned using GMAP (Wu and Watanabe 2005;Wu et al. 2016) to hg19. The GMAP alignments were filtered using quality scores from BWA-MEM (Li 2013) alignments by removing any reads that in BWA-MEM have alignments below a mapping quality of 60. The remaining GMAP alignments were used for gene fusion detection using ToFU (Gordon et al. 2015). Aligned fusion transcripts identified by ToFU were intersected with the GENCODE hg19 annotation (Harrow et al. 2012), and the total number of full-length reads supporting fusions between each pair of genes was counted. All putative gene fusions with at least five full-length Iso-Seq reads from ToFU were input into SplitThreader (Nattestad et al. 2016a) to identify those with any combination of long-range variants that place the genes within 100 kbp of each other. Gene fusion alignments were visualized and figures generated using Ribbon (Nattestad et al. 2016b). The Ribbon source code is also available as a Supplemental File.

Data access
Illumina and PacBio sequencing data from this study have been submitted to the NCBI BioProject (https://www.ncbi.nlm.nih. gov/bioproject) under accession number PRJNA476239. Variant calls from all callers are available as a Supplemental File. The alignments, assemblies, and variant calls are also available at: http:// www.schatz-lab.org/publications/SKBR3/.